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ABSTRACT 


Liquid  cooling  of  a  three-by-three  array  of  commercially  available  leadless  chip  carrier 
packages,  mounted  on  a  ceramic  substrate  was  examined.  Baseline  data  were  obtained  for 
cooling  with  pure  dielectric  liquids.  The  effects  of  addition  of  high  thermal  conductivity  ceramic 
powder  to  the  liquid  were  next  examined,  both  for  natural  and  forced  circulation  conditions. 
Vertical  and  horizontal  orientations  were  studied,  for  two  different  ceramic  particle  types,  and 
two  different  particle  sizes  for  each  ceramic.  For  a  range  of  chip  power  levels,  chip,  substrate  and 
cold  plate  temperatures  were  measured.  Interpretations  for  these  data  are  provided.  A  numerical 
model  was  developed  for  the  vertical  geometry  and  compared  to  the  measurements  obtained. 
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L  INTRODUCTION 


A.  STATEMENT  OF  PROBLEM 

As  the  demand  for  faster  yet  smaller  computers  has  grown,  one  physical  limitation  has 
been  the  heat  removal  capacity  from  the  system.  The  search  for  effective  thermal 
management  techniques  has  expanded  in  recent  years  and  is  expected  to  continue  as  long 
as  package  sizes  decrease  and  clock  speeds  increase.  These  two  factors  combine  to 
dramatically  increase  the  heat  flux  seen  by  the  package. 

Historically,  there  has  been  a  progression  in  the  heat  removal  capacity'  that  has 
matched  or  only  slightly  trailed  the  increase  in  the  heat  generation  rate.  Thus,  thermal 
considerations  have  generally  had  a  secondary  effect  in  slowing  the  quest  for  greater 
computational  speed  and  reduced  package  size.  Today,  though,  the  heat  fluxes  exceed  6S 
watts/cm^  and  the  capability  to  remove  this  heat  flux  is  sorely  challenged.  Bar-Cohen 
[Ref  1] 

Natural  circulation  air  cooling  was  one  of  the  first  methods  used  to  remove  heat.  The 
advantages  are  obvious,  air  is  virtually  firee  and  is  usually  found  in  abundance.  However, 
for  higher  heat  loads  and  for  applications  with  concern  for  airborne  contaminant  damage, 
natural  convection  may  not  be  adequate  or  appropriate.  Forced  convection  air  cooling  can 
provide  increased  heat  removal  capability.  For  many  applications  the  pressure  drop  due  to 
filtering  is  within  acceptable  ranges.  Sloan  [Ref  2] 


Mezenq  et.al.  [Ref.  3]  studied  gas  fluidized  beds  to  improve  heat  transfer  capability 
They  demonstrated  dramatic  performance  increases  (heat  transfer  coefficient  increase)  by 
the  use  of  small  suspended  sand  or  silica  particles  in  an  air  stream.  The  bed  exhibited  a 
very  high  effective  thermal  conductivity  with  low  dependence  on  particle  conductivity 
Situations  arise,  where  either  due  to  the  heat  load,  or  the  environment,  air  is  not  a 
suitable  candidate  as  the  primary  heat  sink.  A  number  of  techniques  are  currently  being 
investigated  for  such  applications.  They  include  conduction  with  indirect  liquid  cooling, 
use  of  fluid  backplane  and  direct  liquid  immersion  cooling. 

One  approach  is  the  use  of  thermally  conductive  solids  to  transfer  the  heat  to  a  cold 
fluid  as  exemplified  by  the  water  cooled  piston  structure  of  the  IBM  Thermal  Conduction 
Module.  A  lower  liquid  side  thermal  resistance  can  be  achieved  by  utilizing  a  highly  finned 
silicon  or  ceramic  heat  exchanger.  The  interfiu;e  resistance  between  the  chip  or  package 
surface  and  the  primary  conductor  poses  the  primary  limitation  on  the  thermal 
performance.  Bergles  and  Bar  Cohen  [Ref  4] 

A  second  approach  is  the  use  of  a  fluid  backplane.  This  can  reduce  the  interfacial 
resistance  between  the  chip  and  the  coolant.  The  resistance  can  be  entirely  eliminated  by 
making  the  cold  plate  an  integral  part  of  the  module.  Kishimoto  and  Osaki  [Ref  S] 

A  third  approach  is  the  use  of  direct  immersion  cooling.  Here,  the  electronic 
components  are  immersed  in  a  dielectric  liquid.  Heat  transfer  in  these  situations  can  be 
either  single  phase  or  with  a  phase  change.  The  coolant  flow  may  be  either  natural  or 
forced.  Various  fluorocarbon  liquids  available  as  Fluorinerts®  (3M  Corporation)  have 
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been  used  in  this  application.  Another  variation  is  the  use  of  mixtures  of  phase  change 
particles  and  dielectric  liquids  to  improve  the  heat  transfer  capabilities  of  the  pure  fluid; 
Choi,  et.al.  [Ref  6]. 


B.  PREVIOUS  RESEARCH  IN  IMMERSION  COOLING 

The  conduction  interfacial  resistance  at  the  chip  or  package  places  an  upper  bound  on 
the  attainable  thermal  performance  of  methods  involving  indirect  liquid  cooling.  Baker 
[Ref  7]  outlined  the  fundamentals  of  immersion  cooling  and  reported  the  results  of  a 
study  using  Freon  -  1 13  and  Dow  Coming  #200  silicone  dielectric  liquid  to  effectively 
cool  small  heat  sources.  Immersion  cooling  is  currently  commercially  implemented  in  the 
Cray-2  supercomputer;  Danielson,  et.  al.  [Ref  8.].  Other  geometries  have  been 
investigated.  One  possible  design  is  a  liquid  encapsulated  module,  similar  to  that  studied 
by  IBM.  In  this  design,  natural  convection  heat  transfer  carries  the  heat  dissipated  from 
the  immersed  chips  to  a  perfluorinated  liquid  and  eventually  to  the  enclosure  walls; 
Bergles  and  Bar  Cohen  [Ref  4].  Related  concepts  involve  pool  boiling  at  the  chips 
investigated  by  Arata  [Ref  9]  and/or  condensation  at  the  enclosure  walls. 

A  number  of  studies  of  natural  convection  in  geometries  of  interest  to  electronic 
cooling  have  recently  been  carried  out.  Joshi  et.al.  [Ref  10]  presented  flow  visualizations 
and  component  surface  temperature  measurements  for  natural  convection  cooling  of  a 
three  by  three  array  of  discrete  heated  protrusions  on  the  vertical  wall  of  a  rectangular 
enclosure  filled  with  a  dielectric  liquid.  Joshi  et.al.  [Ref  1 1]  presented  experimental 
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studies  of  the  heat  transfer  and  flow  characteristics  of  a  column  of  protruding  heat  sources 
on  a  vertical  surflice  and  within  a  vertical  channel.  Sathe  and  Joshi  [Ref  12]  reported 
results  of  a  two-dimensional  numerical  investigation  of  natural  convection  flow  and  heat 
transfer  arising  from  a  protruding  heat  source  on  a  vertical  plate  within  an  enclosure.  Joshi 
and  Paje  [Ref  13]  reported  experimental  results  of  natural  convection  heat  transfer  from  a 
commercially  available  leadless  chip  carrier  package.  Wroblewski  and  Joshi  [Ref  14] 
reported  the  results  of  a  numerical  investigation  of  th>  same  leadless  chip  carrier  package 
studied  by  Joshi  and  Paje. 

The  thermal  properties  of  four  Fluorinerts  and  water  are  shown  in  Table  1-1 . 

Reference  to  Table  1-1  reveals  these  Fluorinerts  to  have  low  thermal  conductivities, 
specific  heats  and  latent  heats  of  vaporization  compared  to  water.  The  direct  liquid 
cooling  of  electronic  components  necessitates  the  use  of  chemically  stable  and  inert, 
non-toxic  liquids  with  high  dielectric  strength  and  high  volumetric  resistivity.  The  need  for 
a  fluid  suitable  for  electronic  circuit  applications  has  constrained  the  choices  to  fluids 
which  have  poor  thermal  transport  characteristics.  The  desire  to  improve  the  transport 
properties  led  to  this  study.  It  examines  the  effect  of  particle  additions  to  the  Fluorinert 
liquids  on  their  heat  transfer  characteristics.  Also  presented  are  measurements  and 
numerical  simulation  of  natural  convection  in  pure  liquids. 
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TABLE  1-1:  THERMOPHYSICAL  PROPERTIES  OF  SEVERAL  LIQUIDS  AT 
ATMOSPHERIC  PRESSURE  (3M  MANUAL)  fREF  IS  ] 


Perfluorinated  Liquid  [>esignation 


Property 

FC-87 

FC-72 

FC-84 

FC-75 

Water 

Boiling  Point, 
degrees  C 

30 

56 

83 

100 

100 

Liquid  Density, 

P,.  kg/m3 

1633 

1680 

1575 

1590 

958 

Kinematic  Viscosity, 

V,  cs 

4.20E-04 

0.4 

0.55 

4.5E-04 

2.70E-04 

Specific  Heat, 
c.  J/kg  K 

1088 

1088 

1130 

1172 

4184 

Thermal  Conductivity, 
k,  W/mK 

5.51E4)2 

5.45E-02 

5.35E-02 

5.70E-02 

6.83E-01 

Vol.  Coef.  Expansion, 

P.  K  ' 

1.60E-03 

1.60E-03 

1.50E-03 

1.40E-03 

2.00E-04 

Dielectric  Constant 

1.71 

1.72 

1.71 

1.75 

78.00 

Average  molecular 
weight,  g/mole 

288 

338 

388 

438 

18 

C.  OBJECTIVES 

The  investigation  reported  here  is  a  continuation  of  the  studies  conducted  by  Joshi  et. 
al.  [Ref  10,  1 1],  Sathe  and  Joshi  [Ref  12],  Joshi  and  Paje  [Ref  13]  and  Wroblewski  and 
Joshi  [Ref  14],  The  present  study  was  conducted  with  a  three  by  three  array  of  leadless 
chip  carrier  packages  mounted  on  a  ceramic  substrate  which  formed  one  wall  of  a 
dielectric  liquid  filled  enclosure.  The  first  part  investigated  the  effect  of  various  ceramic 
particle  loadings  on  overall  heat  transfer  rates  in  two  orientations,  the  (i)  horizontal  and 
(ii)  vertical.  Ceramics,  particularly  the  nitrides,  are  a  class  of  materials  that  combine  many 
useful  features.  They  have  very  high  electrical  resistivities,  moderate  densities  and  have 
exceptionally  high  thermal  conductivities.  Table  1-2  shows  the  properties  of  some 
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ceramics  that  were  considered  for  this  study.  The  present  investigation  studied  the  effect 
of  powdered  Boron  Nitride  (BN)  and  Aluminum  Nitride  (AIN)  additions  on  the  thermal 
transport  characteristics  of  Fluorineil-75  (FC  75).  During  the  course  of  the  investigation, 
no  enhancements  were  noted  in  the  heat  transfer  rates  as  a  result  of  the  particle  additions. 
Extensive  measurements  of  heat  transfer  characteristics  were  made  for  various  particle 
sizes  and  volume  fractions.  The  horizontal  geometry  was  tested  in  natural  circulation, 
natural  circulation  with  external  vibration,  and  forced  circulation  modes.  The  vertical 


geometry  was  only  tested  in  the  natural  circulation  mode. 
TABLE  1-2;  CERAMIC  PROPERTIES 


Property 

Aluminum 

Nitride 

(AIN) 

Boron 

Nitride 

(BN) 

Alumina 

AlA 

Magnesia 

MgO 

Silicon 

Nitride 

(Si,N,) 

Density, 
p,  kg/in3 

3050* 

2290* 

3970* 

3580* 

3440* 

Thermal  Q)iiductivity, 
k,  W/mK 

320* 

225* 

30* 

48* 

16  -  33* 

Electrical  Resistivity 
n-m 

>10E12* 

>10E12* 

>10E12* 

>10E12* 

>10E12* 

1  *  Samsonov  (Ref.  16],  *  Slack  (Ref  17],*Caliister  (Ref  18],  *  Mutterties  (Ref  19] 

The  second  part  of  the  study  was  the  modification  and  use  of  a  numerical  model  to 
simulate  the  three  dimensional  transport  for  the  pure  liquid  conditions  with  a  vertical 
orientation  of  the  substrate. 

The  specific  objectives  were; 

•  To  design  and  build  an  enclosure  to  study  natural  convection  fi'om  a  three  by  three 
package  array  in  fluorocarbon  liquids. 

*  To  study  the  effects  of  component  orientation  on  the  overall  heat  transfer  rates  in 
the  enclosure  under  natural  convection  conditions. 
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To  design  and  build  an  enclosure  in  order  to  study  forced  convection  from  the 
above  package  array  in  fluorocarbon  liquids. 

To  investigate  heat  transfer  in  liquid->ceramic  mixtures  for  several  ceramic  types, 
particle  sizes  and  volume  fractions. 

To  develop  a  numerical  model  of  the  nine  package  vertical  substrate  orientation. 


n.  EXPERIMENTAL  APPARATUS 


AND  PROCEDURE 

As  shown  in  Figure  2-1,  the  experimental  apparatus  consisted  of  several  components; 
an  enclosure  containing  the  electronic  packages,  a  power  distribution  system,  a  heat 
removal  system,  and  a  data  acquisition  system.  For  the  forced  circulation  flow  data,  a 
screw  type  pump,  a  turbine  flow  meter  and  associated  tubing  were  added  to  the  apparatus. 
Figure  2-2  is  a  photograph  of  the  enclosure  interior. 


Figure  2-1:  System  arrangement  and  apparatus. 
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Figure  2-2:  Photograph  of  enclosure  interior. 


A.  APPARATUS  FOR  NATURAL  CONVECTION  TESTS 

The  enclosure  consisted  of  four  side  walls  of  2.S  cm  thick  Plexiglas  plates,  with  an 
outside  dimension  of  9.86  cm  by  9.92  cm,  a  bottom  wall  of  2.5  cm  thick  Plexiglas,  and  an 
aluminum  cover  plate  9.96  cm  by  9.96  cm  and  0.66  cm  thick.  Both  the  enclosure  and  the 
plate  had  a  reference  chamfer  cut  in  the  front  right  comer  to  permit  retention  of  the 
proper  alignment  of  the  enclosure  and  plate.  The  interior  of  the  enclosure  initially 
measured  S.  18  cm  by  S.  12  cm  by  5.08  cm  deep.  Figure  2-3  shows  an  overhead  view  of  the 
interior  of  the  enclosure. 

The  rear  wall  was  provided  with  a  9  pin  connector  to  provide  a  data  and  power  path 
to  the  package  assembly.  A  ceramic  plate  5.08  cm  by  5.08  cm  by  0.05  cm  thick  was 


inserted  to  protect  the  leads  and  to  prevent  flow  disturbances  caused  by  the  wires.  The 
ceramic  plate  reduced  the  interior  to  4.72  cm  by  5.08  cm  by  5.08  cm  deep.  A  recessed 
channel  with  a  3  mm  diameter  gasket  was  provided  in  the  top  of  the  side  walls.  The  rear 
wall  also  had  a  0.64  cm  hole  drilled  to  provide  a  vent/expansion  path.  A  0.64  cm  I  D. 
tygon  tube  was  attached  to  provide  a  surge  volume  during  the  heating  cycle.  Figure  2-4 
shows  the  enclosure  section  view,  cut  through  the  middle  row  of  packages.  Additionally,  a 
form  fitting  2.54  cm  thick  insulation  pad,  which  is  not  shown,  surrounded  the  enclosure 
during  operation. 


Figure  2-3;  Plan  view,  enclosure  intoior,  cover  removed. 


10 


Figure  2-4;  Section  view,  enclosure  midplane. 

Most  of  the  power  generated  by  the  package  was  renmved  from  the  enclosure  by  the 
aluminum  cover  plate.  Eight  0.64  cm  holes  were  drilled  and  tapped  to  provide  a  means  of 
attaching  the  plate  to  the  enclosure.  Four  thermoelectric  (Peltier  effect)  coolers  were 
mounted  to  the  top  of  the  aluminum  plate  to  remove  the  heat  produced  by  the  electronic 
packages.  The  four  coolers  were  in  turn  cooled  by  four  Pin  Grid  Arrays  (PGA)  with 
integral  cooling  fims.  The  four  coolers  were  powered  by  a  HP  6286  DC  Power  Supply 
under  control  of  a  HP  8SB  micro  computo’.  The  coolers  were  controlled  to  maintain  a 
temperature  of  24  °C  on  the  lower  fiu:e  of  the  aluminum  plate.  This  temperature  was 
monitored  1^  four  copper-constantan  thermocouples.  The  Pin  Grid  Array  cooling  fans 
were  also  controlled  by  a  HP  6289A  DC  Power  Supply.  These  Pin  Grid  Arrays  and  their 
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associated  fans  are  normally  installed  in  computers  utilizing  Intel  486  chips  and  were 
operated  in  a  similar  manner. 

The  electronic  packages  formed  a  three  by  three  array  of  leadless  chip  carriers 
mounted  on  a  ceramic  substrate.  The  ceramic  substrate  of  S.  18  cm  length,  5.06  cm  width 
and  a  thickness  of  0.07  cm,  was  mounted  horizontally  at  the  base  of  the  enclosure.  Five  of 
the  nine  packages  were  equipped  with  diode  type  temperature  sensors  integrated  within 
the  chips.  These  diodes  display  a  linear  decrease  in  resistance  with  temperature  increase. 
Only  three  of  the  five  diodes  were  monitored  due  to  space  limitations.  Like  any  solid  state 
component,  the  diodes  are  limited  in  their  range  of  temperature.  The  manu&cturer  (Texas 
Instruments)  indicated  a  maximum  of  ISST,  and  no  data  was  acquired  above  this 
temperature.  The  temperature  sensors  are  mounted  on  the  base  of  the  package.  The  chip  is 
covered  by  a  thin,  brass  lid.  Between  the  chip  and  the  lid  a  small  air  gap  exists.  The  nine 
chips  whlun  the  package  array  could  be  wed  in  a  wide  variety  of  configurations.  Due  to 
space  connderations,  all  nine  chips  were  wired  in  paralld.  The  wire  gage  limited  the  total 
current  to  1  amp  DC.  This  proved  to  be  sufficient  during  these  trials  and  exceeded  the 
heat  removal  capacity  of  the  cooling  medium. 

Attached  to  the  bottom  surfime  of  the  ceramic  substrate,  under  each  package  was  a 
0. 127  mm  diameter  copper-constantan  thermocouple.  These  thermocouples  were 
calibrated  at  the  same  time  as  the  internal  diodes  in  a  temperature  controlled  bath  against  a 
platinum  resistance  thermometer.  Four  0. 127  mm  coppo’-constantan  thermocouples  were 
attached  to  the  bottom  of  the  aluminum  plate.  An  additional  0. 127  mm  copper-constantan 
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thermocouple  monitored  the  temperature  of  the  thermoelectric  coolers  to  prevent  damage 
to  the  coolers.  Two  thermocouples  were  mounted  on  one  enclosure  outside  wall  to  assist 
in  determining  steady  state  conditions.  Ambient  temperature  was  also  monitored. 

B.  STUDIES  WITH  EXTERNAL  EXCITATION 

As  will  be  discussed  later,  the  ceramic  particles,  especially  the  Aluminum  Nitride, 
were  denser  than  the  Fluorinert,  and  were  not  able  to  be  picked  up  by  the  buoyant  stream. 
In  an  effort  to  keep  them  from  settling  on  the  substrate,  an  external  excitation  was 
provided  via  an  eccentric  weight  attached  to  a  DC  motor  via  a  common  bed  plate  to  the 
enclosure.  No  other  modifications  to  the  enclosure  were  made  for  these  data  runs. 

C.  MODfflCATIONS  FOR  FORCED  CONVECTION 

For  the  forced  circulation  portion  of  the  study,  a  flow  loop  was  constructed  to  enable 
fluid  to  be  circulated  through  the  enclosure  as  illustrated  in  Figure  2-S  and  2-6.  Four  0.64 
cm  diameter  holes  were  drilled  through  the  ri^t  wall  of  the  enclosure  used  for  the  natural 
circulation  study.  The  holes  were  arranged  with  two  in  the  lower  portion  of  the  wall  and 
two  in  the  upper  portion.  A  divider  plate  was  inserted  within  the  enclosure  so  as  to  allow 
the  fluid  entering  at  the  lower  two  holes  to  flow  over  the  package,  up  through  three  holes 
drilled  in  the  plate  to  the  uppa*  portion  of  tl^  enclosure,  past  the  cold  plate,  to  the  suction 
of  a  po^ive  displacmient  screw  type  pump,  through  a  turbine  flowmeter  and  return  to  the 
enclosure.  A  hose  clamp  was  applied  to  the  vent  line,  to  prevent  leakage. 
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Figure  2-S;  Plan  view,  enclosure  interior,  cover  removed,  flow  path  illustrated. 


Figure  2-6:  Section  view,  enclosure  midplane,  flow  path  illustrated. 
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D.  POWER  DISTRIBUTION  SYSTEM 

Power  to  the  package  resistors  was  supplied  by  a  0-1.5  A,  0-40  V  HP-6289A  DC 
Power  Supply.  As  discussed  above,  the  individual  packages  were  wired  in  parallel  The 
data  acquisition  system  was  connected  to  read  the  voltage  drop  across  the  nine  packages 
and  also  across  a  known  precision  resistor  that  was  wired  in  series  with  the  package 
assembly.  From  these  two  voltages,  the  input  power  to  the  chipset  using  the  following 
relation; 

Power  =  ^*Vp„k^, 

E.  DATA  ACQUISITION  ASSEMBLY 

The  data  acquisition  assembly  conasted  of  a  HP  300  computer  system,  a  HP-38S2 
Data  Acquisition  Unit  and  a  386  clone  PC.  The  HP-300  computer  instructed  the  data 
acquisition  system  to  monitor  voltages,  resistance  and  temperatures  from  the  desired 
elements.  A  monitoring  program  determined  when  steady  state  had  been  achieved  and 
alerted  the  operator.  The  monitoring  program  also  detected  out  of  limit  parameters, 
similaily  alerting  the  operator  to  take  corrective  action.  When  steady  state,  as  indicated  by 
a  less  than  0.2  °C  change  in  temperatures,  on  the  aluminum  plate,  was  achieved,  a 
separate  data  collection  program  was  run  which  directed  the  data  to  the  386  computer 
where  it  could  be  analyzed.  Following  evalutoion  of  a  successful  data  run,  the  monitoring 
program  was  reloaded  and  the  next  set  of  conditions  established.  For  the  forced 
circulation  data,  an  Omega  Engineering,  Inc.  FLSC-18B  turbine  flowmeter  was  used  to 
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measure  the  flowrate.  The  flowmeter  was  calibrated  using  FC-7S  in  the  range  SI  to  117 
ml/s. 


F.  EXPERIMENTAL  PROCEDURE 

The  study  examined  two  different  geometries,  four  different  fluid/ceramic  mixtures,  as 
many  as  thirteen  different  ceramic  particle  loads  per  ceramic  and  up  to  five  power  levels 
per  ceramic  powder  load,  without  exceeding  the  maximum  chip  temperature  of  135  °C. 

To  prepare  for  each  run,  the  same  basic  procedure  was  followed.  The  enclosure  was  filled 
with  FC  -  7S,  and  the  ceramic  load  was  adjusted  as  necessary.  After  reassembling  the 
cover  and  vent,  the  Pin  Grid  Arrays  were  placed  on  the  top  of  the  Peltier  effect  coolers. 
The  monitoring  system  and  all  auxiliary  power  systems  were  started.  The  monitoring 
program  was  used  during  heatup,  to  check  on  system  performance  and  to  evaluate 
proximity  to  steady  state.  During  the  heatup,  dissolved  air  was  vented  fi’om  the  assembly, 
and  no  data  was  takoi  until  this  was  complete.  A  series  of  runs  was  defined  as  a  set  of 
runs  at  power  levels  decrea^g  fi'om  the  lughest  chosen  for  the  given  conditions;  up  to  a 
maximum  of  five  power  levels  for  each  cenuiuc  load.  The  natural  circulation  and  external 
«(Cttation  runs  all  started  at  1.55  watts  per  chip  (14  watts  total  to  array)  or  the  highest 
achievable  power  levd  and  included  data  at  1.39, 1.22,  l.OS  and  0.89  watts  per  package. 
The  forced  circulation  series  all  started  at  2.2  watts  per  chip  (20  watts  total  to  the  array), 
and  included  data  at  1.55  and  0.9  watts  per  chip. 
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1.  Natural  Circulation  Runs 


Once  the  monitoring  program  indicated  that  steady  state  had  been  achieved,  the 
data  collection  program  (see  Appendix  A)  was  loaded.  No  adjustments  were  pemtitted 
during  data  collection  to  the  input  power  or  thermoelectric  cooler  voltage  Three  data  runs 
were  generally  taken  for  each  condition.  If  the  plate  temperature  or  input  power  varied 
outside  of  the  specified  band  during  the  data  runs,  the  collection  program  was  unloaded, 
and  the  monitoring  program  was  reloaded  to  restore  steady  state  Once  steady  state  was 
reestablished,  data  taking  was  recommenced.  After  satisfactory  runs  were  obtained  for  the 
given  power  level,  the  input  power  was  reduced  to  the  next  point  in  the  series.  The  desired 
power  levels  were  14.0,  12. S,  11.0,  9.S,  and  8.0  watts.  Some  ceramic  loadings  blanketed 
the  packages  and  prevented  achieving  all  the  desired  power  levels. 

2.  Natural  Circulation  -  External  Excitation 

The  runs  with  external  excitation  were  conducted  in  the  same  manner  as  the 
unexcited  natural  circulation  runs,  with  the  exception  that  a  small  direct  current  motor  was 
used  to  vibrate  the  enclosure.  Figure  2-7  is  a  picture  of  the  motor,  eccentric  weight  and 
plate  used  for  the  external  excitation  runs. 

3.  Forced  Circulation 

The  forced  circulation  data  was  collected  generally  in  a  numner  similar  to  the 
natural  circulation  data.  Some  differences  existed.  The  pump  was  used  to  suction  drag 
Fluorinert  and  fill  the  system.  The  small  volume  of  trapped  air  was  vented  through  the 
vent  line.  When  ceramic  was  to  be  added,  the  pump  suction  line  was  disconnected  and  the 


17 


ceramic  was  added  to  the  tubing.  The  pump  was  then  run  to  distribute  the  ceramic  evenly, 
prior  to  energizing  the  chips.  The  ceramic  addition  was  limited  by  the  tendency  of  the 
ceramic  to  clog  the  turbine  flowmeter.  Following  system  fill,  the  monitoring  and  power 
supply  system  were  energized  and  data  was  collected  in  the  same  manner  as  the  natural 
circulation  data,  with  the  same  program.  Following  the  collection  of  a  series  of  data,  a 
separate  program  was  used  to  measure  the  flow  rate,  (see  Appendix  A) 


Figure  2>7;  Motor,  weight  and  plate  used  for  external  excitation  runs. 
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m.  EXPERIMENTAL  RESULTS 


A.  DESCRIPTION  OF  RESULTS 

Contained  in  Appendix  B  is  a  complete  listing  of  all  the  data  runs  performed  for  this 
study.  For  this  investigation,  measurements  were  made  of  chip  temperatures  (using  diode 
sensors ),  ceramic  substrate  temperatures  (using  thermocouples),  cold  plate  temperatures 
(using  thermocouples),  package  voltage  drop  and  precision  resistor  voltage  drop  for 
power  input  determination.  The  data  has  been  collated  by  particle  size,  type  and  heating 
configuration.  For  each  condition,  two  figures  are  plotted,  (i)  total  package  input  power 
versus  the  temperature  difference  between  the  package  diodes  and  the  cold  plate,  (ii)  total 
package  input  power  and  the  peak  diode  temperature.  Figures  3-1  through  3-8  show  the 
results  of  natural  circulation  for  each  of  the  particle  types  in  the  horizontal  geometry. 
Figures  3-9  through  3-12  show  the  results  of  natural  circulation  in  the  vertical  geometry. 
Figures  3-13  and  3-14  show  the  effect  of  external  excitation  for  one  particle  type  and  size. 
Figures  3-1 S  and  3-16  show  the  effect  of  forced  circulation  for  one  particle  type  and  size. 

The  figures  discussed  above  are  in  dimensional  form  and  some  difficulty  arises  when 
trying  to  compare  the  data  from  different  combinations  of  ceramic  and  fluid.  The  results 
were  non-dimensionalized  to  permit  comparison  of  the  data.  During  the  original  testing, 
no  temperature  data  was  taken  directly  on  the  side  of  the  substrate  exposed  to  the  fluid. 
However,  the  temperature  on  the  bottom  of  the  substrate  was  measured  by  nine 
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thermocouples  as  was  the  external  temperature  of  the  case  This  permitted  estimation  of  a 
heat  flux  through  the  enclosure  case,  and  then  a  determination  of  the  substrate  fluid 
interface  temperature.  Once  this  surface  temperature  was  estimated,  a  Rayleigh  number 
was  determined.  The  data  collection  system  provided  sufficient  data  for  the  determination 
of  a  Nusselt  number  directly.  All  of  the  following  figures  are  in  non-dimensional  form. 
Figure  3-17  shows  the  plot  of  Ra  versus  Nu  for  three  different  loadings  of  the  3  micron 
BN  particle  and  a  baseline  run.  Figure  3-18  shows  the  results  of  natural  circulation  for 
each  of  the  different  particle  types  in  the  horizontal  geometry. 

1.  Description  Of  Data  -  Natural  Circulation 

As  previously  discussed,  the  test  surface  consisted  of  a  three  by  three  array  of 
leadless  chip  carrier  packages,  numbered  as  shown  in  Figure  2-3  .  Only  the  center  column 
temperature  sensing  diodes  (labeled  2,  S  and  8  in  Figure  2-3)  were  cormected  to  the  data 
acquisition  system.  All  nine  packages  were  powered  in  parallel,  however,  no  method  was 
available  to  determine  the  variation  in  power  input  between  the  various  packages.  As 
shown  in  Figure  2-4,  nine  thermocouples,  lid)eled  in  the  same  order  as  the  packages  in 
Figure  2-3,  mounted  to  the  base  of  the  ceramic  substrate  were  connected  to  the  data 
acquisition  system.  These  temperatures  and  the  total  input  power  are  represented  in  the 
following  16  figures. 
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Figure  3-2:  Power  versus  package  to  peak  diode  temperature,  horizontal  geometry. 
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Figure  3-4;  Power  versus  padage  to  peak  diode  temperatuie,  horizontal  geometry. 
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Figure  3<5;  Power  versus  package  to  plate  temperature  diflerence,  horizontal  geometry. 


Figure  3-6:  Power  versus  padagc  to  peak  diode  temperature,  horizontal  geometry. 


Figure  3-7:  Power  versus  padcage  to  plate  tenuperature  difference,  horizontal  geometry. 
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Figure  3*9;  Power  versus  package  to  plate  temperature  difiference,  vertical  geometry. 


Figure  3-10:  Power  versos  package  to  peak  diode  temperature,  vertical  geometry. 


Figure  3*11;  Power  versus  package  to  plate  tempenture  difference,  vertical  geometiy. 


Figure  3*12:  Power  versui  package  to  peak  diode  temperature,  vertical  geometry . 
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Figure  3*14:  Rower  versus  purkagr  to  peak  diode  temperature,  external  excitation. 


Power  versus  Temperature  Difference 
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Figure  3>15:  Power  versus  package  to  plate  temperature  difference,  forced  circulation. 


Figure  3-16:  Power  versus  padcage  U)  peak  diode  temperature,  forced  circulation. 
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B.  DISCUSSION  OF  HORIZONTAL  PLATE  NATURAL  CIRCULATION  DATA 


A  baseline  run  was  conducted  using  pure  FC>7S.  This  baseline  was  periodically 
performed  throughout  the  duration  of  data  gathering  to  verify  the  data  taking  system  and 
detect  system  variations.  In  the  horizontal  geometry,  two  different  t>pes  of  ceramic 
particles,  each  having  two  different  sizes  were  investigated.  Boron  Nitride  was  tested  in 
particles  of  0.3-0.7  microns,  hereinafter  referred  to  as  the  0.7  micron  BN  and  2-3  microns, 
hereinafter  referred  to  as  the  3  micron  BN.  Aluminum  Nitride  was  tested  in  particle  sizes 
of  S  microns  and  44  microns.  Table  3-1  shows  the  representative  loadings  that  have  been 

called  'light',  'medium'  and  'heavy'  for  ease  of  reference. 

TABLE  3-1;  PARTICLE  LOADING  NOMENCLATURE 

Horizontal  Geometry 


Light 

^fedillm 

Heavy 

AIN  44.0  micron 

0.1 

1.11 

NA 

AIN  S.O  micron 

1.25 

10.33 

FuU 

BN  3.0  micron 

0.17 

0.39 

Full 

BN  0.7  micron 

O.i 

0.25 

NA 

Vertical  Geometry 


Light 

K/bdium 

Heavy 

AIN  5.0  micron 

5.07 

10.33 

Full 

BN  3.0  micron 

0.17 

1.95 

Full 

External] 

Excitation 

Light 

Medium 

Heavy 

BN  3.0  micron 

0.22 

1.22 

NA 

Forced  Circulation 

Light 

Medium 

Heavy 

BN  3.0  micron 

1.00 

3.66 

NA 

All  weights  in  grams 


1.  Effect  Of  Particle  Loading 

While  many  of  the  particle  types  were  tested  at  various  loadings,  three  general 
classifications,  light,  medium  and  heavy  were  developed  for  the  loading.  Throughout  the 
testing,  the  heavier  loads  were  seen  to  diminish  the  heat  transfer  capability  Figures  3-1 
and  3-2  illustrate  the  point  graphically.  As  the  load  is  increased  to  a  point  where  the 
enclosure  is  packed  as  full  as  possible  with  AIN,  the  peak  power  is  reduced.  At  this  point 
the  input  power  has  been  reduced  to  S0%  of  the  desired  peak  input  power,  yet  the  diode 
temperature  has  risen  by  30°C  above  the  baseline  temperature  at  14  watts  total  input 
power. 

2.  Effect  Of  Particle  Siie 

As  was  discussed  earlier,  two  differtmt  particle  sizes  were  tested  for  each  of  the 
two  ceramics.  Two  pairs  of  figures,  Figures  3-1  and  3-3  and  Figures  3-5  and  3-7  show  the 
effect  of  two  different  particle  sizes.  The  first  figure  of  each  pair  (Figure  3-1  and  3-S) 
represents  the  smaller  of  the  two  particles  for  each  material,  while  the  second  figure  of 
each  pair  (Figure  3-3  and  3-7)  represents  the  larger  of  the  two  particles  for  each  material. 
While  the  loadings  are  not  exactly  the  same,  the  larger  particles  show  much  more  rapid 
heat  transfisT  decranent  as  the  loading  is  increased.  This  was  further  confirmed  by  visual 
observation,  as  the  larger  particles  were  not  picked  up  by  the  buoyant  streams  at  lower 
loadings.  Figures  3-2,  3-4, 3-6  and  3-8  show  the  same  trend  even  more  prominently  in  the 
rapid  rise  of  peak  diode  temperature  for  given  power  levels.  The  larger  particles  were  so 
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deleterious  that  the  desired  input  power  was  not  able  to  be  obtained.  For  these  conditions 
the  input  power  was  limited  to  avoid  exceeding  a  diode  temperature  of  13S,°C 

3.  Effect  Of  Particle  Density 

Boron  Nitride  has  a  density  relatively  close  to  that  of  FC-7S  and  was  observed  to 
be  easily  lifted  off  the  plate  by  the  plumes  that  formed  above  the  chips.  Aluminum  Nitride 
is  quite  a  bit  denser  than  Boron  Nitride  and  as  loading  was  increased,  rapidly  blanketed  the 
plate  and  could  not  be  lifted  by  the  buoyant  plumes.  Figures  3-1  and  Figures  3-5  show  that 
as  the  Aluminum  Nitride  load  is  increased,  the  ability  to  transfer  heat  is  reduced  faster  than 
the  capability  is  reduced  for  the  Boron  Nitride.  In  fact  for  the  heavy  load  with  Boron 
Nitride,  the  peak  package  temperature  is  nearly  ten  degrees  cooler  than  for  the  Aluminum 
Nitride  case.  This  is  quite  unexpected  as  bulk  Aluminum  Nitride  has  a  thermal 
conductivity  nearly  twice  that  of  Boron  Nitride. 

C.  DISCUSSION  OF  VERTICAL  SUBSTRATE  ORIENTATION  NATURAL 
CIRCULATION  DATA 

A  baseline  set  of  data  was  obtained  with  the  enclosure  mounted  with  the  substrate 
containing  the  packages  in  a  vertical  arrangement.  The  cold  plate  was  thus  also  mounted 
vertically,  opposite  the  packages.  This  geonwtry  was  only  tested  in  a  limited  number  of 
runs  as  it  was  observed  that  the  particles  soon  fdl  out  of  suspension  and  rested  on  the 
bottom  wall. 
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1.  Effect  Of  Particle  Loading 

The  Boron  Nitride  3  micron  particles  and  the  Aluminum  Nitride  5  micron  particles 
were  the  only  types  tested  across  a  wide  range  of  particle  loads.  Figures  3-11  and  3-12 
illustrate  that  the  particles  had  a  tendency  to  fall  out  of  suspension.  The  temperatures  soon 
came  to  the  same  values  as  the  baseline  data  when  there  were  insufficient  particles  to  rest 
against  the  vertical  plate.  Figures  3-9,  3-10,  3-1 1  and  3-12  demonstrate  that  when  there 
were  suflk  ient  particles  in  contact  with  the  plate,  reduction  in  the  heat  transfer  occurred. 

2.  Effect  Of  Geometry 

Comparison  of  Ftgu;e>  3-1  and  3-2  with  3-9  and  3-10,  demonstrates  that  the 
vertical  geometry  had  a  better  heat  transfer  coefficient  than  the  horizontal  geometry.  This 
is  to  be  expected,  since  the  horizontal  geometry  develops  cells  that  interfere  with  their 
neighbors.  The  vertical  geometry  results  in  a  larger  circulation  that  sweeps  the  entire  wall. 

D.  DISCUSSION  OF  EXTERNALLY  EXCITED  NATURAL  CIRCULATION 
DATA 

During  the  early  parts  of  the  study  it  was  noted  that  as  the  particle  loading  increased, 
the  fraction  of  the  particles  that  were  lifted  the  buoyant  plumes  decreased  and  the 
packages  became  covered  by  particles,  decreasing  the  heat  transfer  and  raising  their 
temperatures.  It  was  also  noted  earlier  that  Boron  Nitride  has  a  density  close  to  that  of  the 
Fluorinert.  An  external  excitation  was  provided  to  the  enclosure  by  a  small  direct  current 
motor  to  keep  the  particles  suspended  in  the  liquid.  As  Figures  3-13  and  3-14  illustrate. 
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there  was  still  a  decrement  in  performance.  The  decrement  was  greater  than  the  stationary 
performance  for  the  same  load. 


E.  DISCUSSION  OF  FORCED  CIRCULATION  RESULTS 

A  baseline  run  was  conducted  using  pure  FC-75.  In  the  forced  circulation  geometry, 
only  the  3  micron  Boron  Nitride  was  investigated.  For  the  forced  circulation  portion  of  the 
study,  flowmeter  voltage,  (converted  to  flowrate)  was  also  measured.  The  apparatus  failed 
after  testing  was  completed  for  one  Reynolds  number. 

As  illustrated  in  Figures  3-15  and  3-16,  the  addition  of  particles  did  not  enhance  the 
heat  transfer  in  the  forced  circulation  condition.  The  different  particle  loadings  exhibited 
essentially  the  same  performance.  The  measurements  were  limited  by  the  ability  of  the 
turbine  flowmeter  to  pass  the  heavier  particle  loadings  and  as  noted  above,  the  apparatus 
failed  due  to  mechanical  cracking  of  the  ceramic  insert  at  the  rear  wall  early  in  this  testing 
regime. 

F.  NONDIMENSIONAL  GRAPHS 

The  data  was  non-dimensionalized  as  discussed  earlier,  with  some  modifications. 
Additional  thermocouples  were  installed  after  the  completion  of  the  originally  scheduled 
testing.  TlMse  were  attached  to  the  tops  of  packages  1, 2, 4  and  5  and  were  used  to  obtain 
an  estinute  of  the  fluid-surftice  interface  temperatures.  At  low  loadings,  the  surface 
temperature  of  the  packages  is  within  several  degrees  of  the  substrate  bottom  temperature. 
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At  a  higher  particle  loading  the  difference  approaches  10°C  at  14  watts  total  input.  The 
data  did  not  meet  the  strict  repeatability  of  the  previous  data  due  to  the  cracking  noted 
above  and  was  used  to  provide  an  estimate  of  adjustment  of  the  substrate  bottom 
temperature  in  conjunction  with  the  technique  described  earlier. 

Only  the  natural  circulation  data  is  presented  in  non  •  dimensional  form,  as  the 
other  modes  show  similar  trends.  Figure  3-17  represents  the  data  from  only  the  3  micron 
Boron  Nitride  data  runs.  Figure  3-18  is  a  representation  of  the  data  from  representative 
samples  for  all  the  ceramic  particles.  It  is  evident  that  for  all  particles  as  the  loading  is 
increased  the  non  dimensional  Nusselt  number  decreases  for  a  given  Rayleigh  number.  As 
the  loading  increased,  the  heat  transfer  capacity  of  the  FC-7S  diminished  substantially. 
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Figure  3>17;  Ra  versus  Nu  for  3  micron  BN 


Figure  3-18:  Ra  versus  Nu  for  natural  circulation  runs. 
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IV.  NUMERICAL  MODEL 


A  numerical  model  of  the  natural  convection  in  the  enclosure  for  a  three  by  three  array 
of  packages  in  the  vertical  substrate  orientation  was  developed  utilizing  a  control  volume 
approach  as  described  by  Patankar  [Ref  20].  The  model  is  similar  to  that  discussed  by 
Wroblewski  and  Joshi  [Ref  14]  for  a  single  package  and  included  the  following  features; 
control  volumes  for  velocities  that  are  staggered  with  respect  to  those  for  temperature  and 
pressure;  a  power  law  scheme;  harmonic  mean  formulation  for  the  interface  diflusivities; 
the  SIMPLER  algorithm  for  the  velocity  pressure  coupling;  and  a  fully  implicit  forward 
difference  scheme  in  time.  The  conjugate  conduction  in  the  solid  domains  was  handled 
numerically  by  solving  the  same  full  set  of  momentum  and  energy  equations  throughout 
the  entire  enclosure,  but  with  a  large  value  of  viscosity  specified  for  the  solid  regions. 


A.  GOVERNING  EQUATIONS 

The  non-dimensional  governing  equations  for  the  three  dimensional  unsteady 
problem,  assuming  constant  properties  and  the  Boussinesq  approximation,  are  as  follows; 
X  momentum; 

4.  a(vu)  a(wu)  _  _«p  . /Ptxw  *  /«iu  .  ,  ajux 

y  momentum; 

av  .  a(uv)  .  acw)  .  a(wv)  _  ./Prxi/i^/aJv  .  a>v  .  a»Vx  .  « 
ar  +  “^  +  ~  +  “« - av+Viu/ 
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z  momentum; 

aw  ,  agjw)  .  a(vw)  acww)  ap  rPrvi/? 
dt  ax  av  az”  az 


^a^w  .  a^w  ,  a^Wv 


energy  (fluid); 

» 4.  £[H®>  4.  £lX2i  +  « /£i  4.  iii  4.  iii'k 

at  ax  av  az  “VprR,;  W  av^  az^' 


energy  (chip); 

r  —  =  c— 1— « fD  .  iii  .  .  1  1 

®at  VprR**  ^  av^  ^az^^^  Hctr 

energy  (unheated  regions) 


r  —  -(  *  *  TR Y— +  ~  +  — ^1 

at  vmu)  W  ^  av’  ^ 


The  appropriate  non-dimensional  parameters  are  Ra  =  gPQlVAavk^  Pr  =  v/a,  U  = 


uAJ^  V  =  v/U„,  W  =  w/U„.  U„  =  (gpQ/k,)''^.  t  =  tU„/l,  0  =  (T-T,)/(Q/Uc,).  P  =  p/pUo^ 


X  = 


x/l,  Y  =  x/1,  Z  =  z/1,  H.=  h^,  4=  lyi,  R,=  Ri=  Ic/k^  (pCp)/(pCp)f  and  C,= 
(P<:p)AP<^p)r  energy  equation  for  unhemed  regions  is  applicable  to  all  of  the  regions 

within  the  package  except  the  chip  itself  The  subscript  i  refers  to  these  various  regions;  i 
=  s  to  the  ceramic  substrate;  i  =  p  to  the  ceramic  package  itself,  i  =  1  to  the  lid,  i  =  g  to  the 
solder,  i  =  r  to  the  air  gap  between  the  chip  and  the  lid,  and  i  =  m  to  the  gold  and  tungsten 
coating. 

The  boundary  conditions  for  the  enclosure  walls  are  as  follows; 

X  =  O;50/aX  =  O,U»O,  V  =  0,W*0 

X  =  Xl;0*O,U  =  O,V  =  O,W  =  O 

Y*  0,  X^;  a0/5Y  -  0,  U  =  0.  V  =  0,  W  =  0 

z«  0,  Xl;  a0/az  »  0,  u = 0,  V = 0.  w = 0 

vdiere  X/1. 
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Figure  4-1  illustrates  the  axis  selection  used  for  the  analysis.  The  substrate  is 
vertically  aligned.  The  axis  normal  to  the  package  is  designated  the  X  -  axis,  the  axis  in  the 
vertical  plane  designated  the  Y  -  axis,  and  the  axis  parallel  to  the  substrate  surface  is 
designated  as  the  Z  -axis. 


B.  SOLUTION  TECHNIQUE 

The  solution  was  obtained  throughout  the  entire  enclosure.  The  initial  nondimensional 
time  step  was  selected  to  be  At  =  1  for  Ra  =  1 . 1 SE3.  An  algorithm  within  the  program 
expanded  the  time  steps  if  convergence  was  reached  within  a  time  step  in  three  iterations. 
Two  iterations  per  time  step  is  the  minimum  required  due  to  the  technique  used  for 
convo-gence  checking. 
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1.  Results 


a.  Base  Case 

The  initial  case  was  similar  to  the  geometry  and  conditions  developed  by 
Wroblewski  and  Joshi  [Ref.  14]  for  the  single  chip  ieadless  package.  The  solver  used  for 
that  study  has  been  ported  from  an  Amdahl  mainframe  and  refined  to  run  on  a  Sun 
Microsystems  Sparc  10  workstation. 

b.  Current  Case 

Once  the  solver  had  been  verified,  the  geometry  was  revised  to  reflect  the  three 
by  three  array  of  Ieadless  chip  carrier  packages  that  were  used  for  the  experimental  work. 
A  low  Rayleigh  number  (Ra  ==1.1  SE3)  was  used  to  test  the  stability  of  the  numerical 
model.  For  FC-7S,  this  Rayleigh  number  is  equivaloit  to  an  input  power  of  only  4.4 
microwatts.  Thus  the  problem  is  primarily  conduction  heat  transfer  and  the  fluid  develops 
very  low  velocities.  Figure  4-2  shows  the  isotherms  for  the  nine  package  geometry  in  the 
Z  -  Y  plane  passing  through  the  plane  adjacoit  to  the  package  lid  surfaces,  at  r  =  959. 
Figure  4-3  is  a  plot  of  the  fluid  velocities  in  the  same  plane  as  the  temperature  plot  of 
Figure  4-2.  It  shows  the  fluid  developing  upward  velocities  in  the  heated  regions.  This  is 
still  in  the  transient  period  discussed  by  Wroblewski  and  Joshi.  The  temperature  contours 
for  a  vertical  cut  in  the  X  -  Y  plane  through  the  center  of  the  middle  column  of  packages 
are  shown  in  Figure  4-4.  The  plumes  are  starting  to  develop  above  the  packages,  but  there 
is  very  little  flow  developing  in  the  enclosure.  Figure  4-S  shows  the  velocities  in  the  same 
plane  as  Figure  4-4. 
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Ra  =  1  15E3 


Figure  4-2;  Isothenns  for  Ra  »  1. 15E3.  in  the  Y  -  Z  plane  at  X  ~  0.29.  t  =  9S9. 


Ra-1 1SE3 


Figure  4-3:  ZY  velocity  vectns  for  Ra  >  1.1SE3,  in  the  Y  -  Z  plane  at  X  «  0.29.  t  =  9S9. 
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The  model  was  now  reprogrammed  for  the  case  of  a  much  higher  Rayleigh 
number.  The  highest  Rayleigh  number  attained  in  this  study  was  1 . 15E7  This  is  equivalent 
to  a  power  input  of  44  milliwatts.  The  following  figures  are  taken  at  the  same  positions  as 
the  corresponding  figures  for  the  case  of  Ra  =  1. 15E3  Due  to  the  temperature  scaling,  the 
nondimensional  temperatures  are  lower,  while  the  actual  temperatures  are  higher  Figure 
4-6  shows  the  isotherms  for  the  nine  package  geometry  in  the  Z  -  Y  plane  passing  through 
the  plane  adjacent  to  the  package  lid  surfaces,  at  t  =  3619.  Figure  4-7  is  a  plot  of  the  fluid 
velocities  in  the  same  plane  as  the  temperature  plot  of  Figure  4-6.  It  shows  the  fluid 
developing  upward  velocities  in  the  heated  regions.  The  temperature  profiles  for  a  venical 
cut  in  the  X  -  Y  plane  through  the  center  of  the  middle  column  of  packages  is  shown  in 
Figure  4-8.  Figure  4-9  shows  the  velocities  in  the  same  plane  as  Figure  4-8.  The  plumes 
are  starting  to  develop  above  the  packages,  and  the  flow  velocities  are  larger  in  the 
enclosure.  In  addition  the  return  flow  from  the  cold  plate  has  begun  and  flow  is  beginning 
to  develop  the  circulation  expected  for  this  type  of  enclosure. 
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Ra  =  1  1SE7 


Figure  4'6;  Isotherms  for  Ra  =  1. 15E7,  in  the  Y  •  Z  plane  at  X  =  0.29.  x  =  3S 16. 


R8  =  1  15E7 


Figure  4*7:  ZY  vdodty  vectors  for  Ra  *  1.1SE7,  in  the  Y  -  Z  ptane  at  X  «  0.29.  x  =  3S16 
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Figure  4-8;  Isotherms  for  Ra  =  1. 15E7  in  the  X-Y  plane  at  Z  =  2.55.  t  =  3516. 


R8  =  1  15^7 
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Figure  4-9:  XY  vdocity  vectors  for  Ra  *  1.I5E7  in  the  X-Y  plane  at  Z  »  2.55.  t  =  3516. 
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V.  CONCLUSIONS 


The  desire  to  improve  the  heat  transfer  performance  of  FIuorinert-75  through  the  use 
of  ceramic  particles  was  unsuccessful  The  particles  and  loadings  tested  all  resulted  in 
unchanged  or  reduced  heut  transfer  capability  than  the  pure  FC-7S.  The  larger  particles 
caused  a  greater  decrement  than  the  smaller  particles  due  to  their  increased  propensity  to 
resist  the  buoyant  streams  and  rest  on  the  face  of  the  chips.  The  3  micron  Boron  Nitride, 
which  was  the  smallest  particle  and  has  the  density  closest  to  the  Fluonnert,  caused  the 
smallest  decrement  in  heat  transfer. 

The  external  excitation  proved  to  be  successful  in  suspending  the  particles  from  the 
plate,  but  still  resulted  in  a  heat  transfer  decrement. 

The  forced  circulation  data  showed  the  relative  insensitivity  of  the  heat  transfer 
capability  of  the  Fluorinert  to  the  addition  of  particles  in  this  regime.  It  also  demonstrated 
the  dangers  of  using  forced  circulation  in  direct  electronic  cooling.  Throughout  the  natural 
convection  testing,  extending  over  a  period  of  greater  than  six  months,  no  damage  had 
beer  encountered  to  the  packages.  In  the  forced  circulation  testing  the  packages  failed  in 
four  days,  most  probably  due  to  flow  induced  vibration  of  the  divider  plate.  The  danger  in 
actual  electronic  equipment  of  similar  phenomena  should  be  a  key  consideration  in  the 
design  of  direct  immersion  cooling  schemes.  The  results  of  this  testing  indicate  that  for  the 
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materials  tested,  and  possibly  for  a  wider  range  of  ceramic  materials,  the  particles  do  not 
enhance  the  heat  transfer  characteristics. 

The  three  dimensional  numerical  model  was  successfully  modified  to  model  the 
geometry  encountered  in  this  study  for  the  vertical  substrate  orientation.  The  model  was 
able  to  adequately  represent  the  expected  temperature  profiles. 
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VL  RECOMMENDATIONS 


The  following  recommendation  are  made  for  further  study: 

1.  Investigate  the  physical  reason  for  the  failure  of  high  thermal  conductivity  particles 
to  raise  the  overall  mixture  thermal  conductivity. 

2.  Investigate  particles  that  more  closely  match  the  density  of  the  Fluorinert. 

3.  Refine  the  mesh  used  for  the  numerical  model  to  concentrate  in  the  regions  with 
large  gradients. 
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APPENDIX  A 

POWER  AND  TEMPERATURE 
ACQUISITION  PROGRAMS 
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A.  DATA  MONITORING  PROGRAM 


100! 

110! 

120! 

130 

140 

150 

160 

170 

180 

190 

200 

210 

220 

240 

250 

260 

261 

263 

270 

280 

290 

300 

310 

320 

330 

340 

350 

360 

370 

380 

390 

400 

410 

420 

430 

440 

450 

460 

470 

480 

490 


!!!!!!!!!!!!!!!!!!! 

DATA  COLLECTION 
MONITORING  FOR  EQUILIBRIUM 
V(0)  MONITORS  VOLTAGE  DROP  ACROSS 
THE  PRECISION  RESISTOR 

V(l)  MONITORS  THE  VOLTAGE  DROP  ACROSS  THE 
CHIP  SET 

0()  MONITORS  THE  DIODE  OHMS 
T()  ARE  THE  TEMPERATURES 
T(0)  IS  THE  AMBIENT 
T(l-9)  ARE  THE  CHIPS 
T(10)  IS  NOT  USED 
T(1 1-14)  ARE  THE  UNDER  PLATE 
T(16-17)  ARE  THE  SIDE  WALLS 
REALV(1),0(2).T(17) 

REAL  Vt(  1 , 1 0),Ot(2, 1 0),Tt(  1 7. 1 0) 

PRINT  "PWR  * 

INPUT  Pwr 
PRINTER  IS  CRT 
CLEAR  SCREEN 
J=0 

OUTPUT  709;''RST" 

OUTPUT  709;''REAL  V(1),0(2),T(16)’' 

OUTPUT  709;"CONFMEAS  DCV,203-204,USE  700" 
ENTER  709;  V(*) 

PRINT  V(0),V(1),V(0)*V(1)/.101 
Vt(0,J)=V(0) 

Vt(l,J)=V(l) 

Vt(0,10)*Vt(0,10)+V(0) 

Vt(l,10)=Vt(l,10>+V(l) 

OUTPUT  709;"CONFMEAS  OHM,200-202,USE  700" 
ENTER  709;0(*) 

0t(0,j)=(7234.555-0(0))/19.9904 
0t(l.j)=(7241.029-0(l))/20.0697 
Ot(2,J)=(724 1 .567.0(2))/20.087I 
PRINT  Ot(0,J) 

PRINT  Ot(l,J) 

PRINT  Ot(2,J) 

Ot(0, 10)==Ot(0, 10)+Ot(0,J) 

Ot(l,10)*Ot(l,10)+Ot(l,J) 

0t(2,10)«0t(2,10)+0t(2,j) 
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500  OUTPUT  709;"CONFMEAS  TEMPT,300-309.3 1 1-3 16,3 18, USE  700  INTO  T" 
510  OUTPUT  709;"VREAD  T" 

520  FOR  1=0  TO  16 
530  ENTER  709;T(I) 

540  IF  I<  1 0  OR  I>  14  THEN  560 
550  IF  T(I)>60  THEN  BEEP 
560  PRINT  I,T(I) 

570  Tt(I,J)=T(I) 

580  Tt(I,10)=Tt(I,10)+T(I) 

590  NEXT  I 
600  Tave=0 
610  FOR  1=10  TO  13 
620  Tave=Tave+T(I) 

630  NEXT  I 
640  PRINT  Tave/4 
650  J=J+1 

660  IF  J=  10  THEN  680 
670  GOTO  320 
680  PRINT  Vt(0,10)/J 

690  PRINT  Vt(  1,1 0)/J 

69 1  Power=Vt(0, 1 0)*  Vt(  1 , 1 0)/(.  1 0 1  *J*J) 

700  PRINT  Power 

710  Vt(0.I0)=0 
720  Vt(1.10)=0 
730  PRINT  Ot(0,I0)/J 
740  PRINT  Ot(l,10)/J 
750  PRINT  Ot(2,10)/J 
760  Ot(0,l0)=0 
770  Ot(l,10)=0 
780  Ot(2,10)=0 
790  FOR  1=0  TO  16 
800  PRINT  I,Tt(I,J)/10 
810  Tt(I,J)=0 

820  NEXT  I 

830  PRINT  "Tave=";Tave/4 

840  PRINT 

84 1  IF  Tave<23 .8*4  OR  Tave>24.2*4  THEN  850 

842  IF  ABS(Power-Pwr)>.05  THEN  850 

843  BEEP 

844  PRINT  "STEADY  STATE  MET" 

850  K=K+1 

860  PRINT  K 
870  PRINT 
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880  GOTO  290 
890  END 


B.  DATA  ACQUISITION  PROGRAM 

100!!!!!!!!!!!!!!!!!!!!! 

1 10!!  DATA  COLLECTION  FILE 

120!!  ENTER  THE  RUN  TYPE 

150  !!  V(0)  MONITORS  VOLTAGE  DROP  ACROSS 

160  !!  THE  PRECISION  RESISTOR 

170  !!  V(I)  MONITORS  THE  VOLTAGE  DROP  ACROSS  THE 
180  !!  PACKAGES 

190  !!  0()  MONITORS  THE  DIODE  OHMS 

200  ! !  TO  ARE  THE  TEMPERATURES 

210  !!  T(0)  IS  THE  AMBIENT 

220  !!  T(l-9)  ARE  THE  CHIPS 

230  !!  T(10)  IS  NOT  USED 

240  !!  T(1 1-14)  ARE  THE  UNDER  PLATE 

250  !!  T(15)  IS  THE  PELTIER 

260  !!  T(16-17)  ARE  THE  SIDE  WALLS 

270  REAL  V(l),0(2),T(17),Ts 

280  REALVt(l.l0),Ot(2,10),Tt(17,10) 

290  PRINTER  IS  CRT 
300  CLEAR  SCREEN 

3 10  PRINT  "INPUT  TARGET  POWER  LEVEL: " 

320  INPUT  Pwr 
330  Geom$="H  -  FC" 

340  PRINTER  IS  9 
350  Fluid$="FC  75" 

360  Ceramic$»"3m  BN" 

370  Grains=3.66 

380  PRINT  "THIS  RUN  IS  AT  ";Pwr,"  WATTS" 

390  PRINT  DATE$(TIMEDATE),TIME$(TIMEDATE) 

400  PRINT  "THE  FLUID  IS:  "^HuidS;",  ARRANGEMENT  IS:  ";Geom$ 
410  PRINT  "THE  CERAMIC  IS:  ";Ceraimc$ 

420  PRINT  "THE  CERAMIC  WEIGHT  IS:  ";Grains;"  GRAMS" 

430  J=0 

440  OUTPUT  709;"RST" 

450  OUTPUT  709;"REAL  V(1),0(2)J(16)" 


460  OUTPUT  709;XONFMEAS  DCV.203-204.USE  700" 

470  ENTER  709; V(*) 

480  PRINT  V(0) 

490  PRINT  V(l) 

500  PRINT  V(0)*V(1)/.  101 

5 10  OUTPUT  709;"CONFMEAS  OHM,200-202,USE  700" 

520  ENTER  709;O(*) 

530  PRINT  (7234  555-O(0))/19.9904 
540  PRINT  (724 1 .029-O(  1  ))/20 .0697 
550  PRINT(7241.567-O(2))/20.087I 

560  OUTPUT  709;''CONFMEAS  TEMPT,300-309.3 1 1-3 16,3 18,USE  700, INTO  T" 
570  OUTPUT  709;"VREAD  T" 

580  FOR  1=0  TO  16 
590  ENTER  709;T(I) 

600  IFI<10THEN640 
610  IF  T(I)>60  THEN  BEEP 
620  IF  I>13  THEN  640 
630  Ts=Ts+T(I) 

640  PRINT  T(I) 

650  NEXT  I 
660  J=J+1 

670  IFJ=10THEN690 
680  GOTO  460 

690  PRINT  DATE$(TIMEDATE),'nME$(TIMEDATE) 

700  PRINT  "DATA  RUN  ENDED" 

710  PRINTER  IS  CRT 
720  PRINT  Ts/40 
730  END 


C.  FLOWRATE  COLLECTION  PROGRAM 

100!!!!!!!!!!!!!!!!!!!!! 

1 10!!  FLOW  METER  CALIBRATION 

120  !!  F(0)  MONITORS  THE  VOLTAGE  ACROSS  THE 

130 !!  FLOWMETER 

140  REALF(0) 

150  REALFt(10) 

160  PRINTER  IS  CRT 
170  CLEAR  SCREEN 
180  PRINT  "RHEOSTAT  SETTING" 

190  INPUT  R 


52 


200  J-0 

210  OUTPUT  709;"RST" 

220  OUTPUT  709;"REAL  F(0)" 

230  OUTPUT  709;"CONFMEAS  DCV,205,USE  700" 
240  ENTER  709;F(*) 

250  F(0)=F(0)*(-1  0) 

260  PRINT  "FLOW  ";F(0) 

270  Ft(J)=F(0) 

280  FOR  N=1  TO  2000 
290  NEXTN 
300  J=J+l 

310  IFJ=10THEN330 

320  GOTO  210 

330  PRINT  "ELAPSED  TIME" 

340  INPUT  Ti 

350  PRINT  "ML  PUMPED" 

360  INPUT  Ml 

370  PRINTER  IS  9 

380  PRINT  "CALIBRATION" 

390  PRINT  "DIAL  SETTING  ".R 
400  FORJ=lTO10 
410  PRINT  Ft(D 
420  NEXTJ 

430  PRINT  "ELAPSED  TIME  ";Ti 
440  PRINT  "VOLUME  (ML)  ";M1 
450  PRINTER  IS  CRT 
460  END 
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APPENDIX  B 
LIST  OF  RUNS 
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Run 

Name 


C0_12.5 


C0_11 


C0_9.5 


Ceramic 


None 


None 


None 


None 


None 


AIN 


AIN 


AIN 


AIN 


AIN 


AIN 


AIN 


AIN 


AIN 


AIN 


AIN 


AIN 


AIN 


AIN 


AIN 


AIN 


AIN 


AIN 


AIN 


B1_9.5 


B1_8 


Size 

Microns 

Load 

Grams 

Power 

Watts 

Orientation 

NA 

NA 

14 

Horizontal 

NA 

NA 

12.5 

Horizontal 

NA 

NA 

11 

Horizontal 

NA 

NA 

9.5 

Horizontal 

NA 

NA 

8 

Horizontal 

5 

1.25 

14 

Horizontal 

5 

1.25 

12.5 

Horizontal 

5 

1.25 

11 

Horizontal 

5 

1.25 

9.5 

Horizontal 

5 

1.25 

8 

Horizontal 

5 

2.52 

14 

Horizontal 

5 

2.52 

12.5 

Horizontal 

5 

2.52 

11 

Horizontal 

5 

2.52 

9.5 

Horizontal 

5 

2.52 

8 

Horizontal 

5 

5.07 

11.5 

Horizontal 

5 

5.07 

11 

Horizontal 

5 

5.07 

9.5 

Horizontal 

5 

5.07 

8 

Horizontal 

5 

10.33 

11.3 

Horizontal 

5 

10.33 

11 

Horizontal 

5 

10.33 

9.5 

Horizontal 

5 

10.33 

8 

Horizontal 

5 

fuN 

8 

Horizontal 

1.95 

14 

Horizontal 

3 

1.95 

12.5 

Horizontal 

3 

1.95 

11 

Horkontal 

1.95 

9.5 

Horizontal 

3 

1.95 

8 

Horizontal 

Run 

Name 

Ceramic 

Size 

Microns 

Load 

Grams 

Power 

Watts 

Orientation 

B2_14 

BN 

3 

.166 

14 

Horizontal 

BN 

3 

.166 

12.5 

Horizontal 

B2_11 

BN 

3 

.166 

11 

Horizontal 

B2_9.5 

BN 

3 

.166 

9.5 

Horizontal 

B2_8 

BN 

3 

.166 

8 

Horizontal 

B3_14 

BN 

3 

1 

0.303 

14 

Horizontal 

B3_125 

BN 

3 

0.303 

12.5 

Horizontal 

B3_11 

BN 

3 

0.303 

11 

Horizontal 

B3_9.5 

BN 

3 

0.303 

9.5 

Horizontal 

B3_8 

BN 

3 

0.303 

8 

Horizontal 

B4_14 

BN 

3 

0.393 

14 

Horizontal 

B4_125 

BN 

3 

0.393 

12.5 

Horizontal 

B4_11 

BN 

3 

0.393 

11 

Horizontal 

B4_9.5 

BN 

3 

0.393 

9.5 

Horizontal 

B4_8 

BN 

3 

0.393 

8 

Horizontal 

BN 

3 

0.016 

14 

Horizontal 

B5_125 

BN 

3 

0.016 

12.5 

Horizontal 

B5_11 

BN 

3 

0.016 

11 

Horizontal 

B5_9.5 

BN 

3 

0.016 

9.5 

Horizontal 

B5_8 

BN 

3 

0.016 

8 

Horizontal 

B6_14 

3 

0.054 

14 

Horizontal 

B6_125 

BN 

3 

0.054 

12.5 

Horizontal 

BN 

3 

0.054 

11 

Horizontal 

B6_9.5 

mm 

3 

0.054 

9.5 

Horizontal 

B6_8 

BN 

3 

0.054 

8 

Horizontal 

B7_14 

■m 

3 

0.097 

14 

Horizontal 

B7_125 

BN 

3 

0.097 

12.5 

Horizontal 

B7_11 

BN 

3 

0.097 

11 

Horizontal 

B7_9.5 

3 

0.097 

9.5 

Horizontal 

B7_8 

BN 

3 

0.097 

8 

Horizontal 

Run 

Ceramic 

Name 

B8_14 


B8_125 


B8_11 


B8_9.5 


D1_14 


D1_125 


D1_11 


D1_9.5 


D2_14 


D2_125 


D2_11 


D2_9.5 


D2_8 


D3_14 


D3_125 


D3_11 


D3_9.5 


E1_14 


E1_125 


E1_11 


E1_9.5 


E1_8 


0.115 


0.115 


0.115 


0.115 


0.115 


BN 

3 

BN 

3 

BN 

3 

BN 

3 

BN 

BN 

■■i 

BN 

3 

BN 

3 

BN 

3 

BN 

BN 

3 

BN 

3 

mm 

3 

BN 

3 

BN 

3 

BN 

3 

BN 

0.7 

BN 

0.7 

BN 

0.7 

BN 

0.7 

BN 

0.7 

0.099 


0.099 


0.099 


0.155 


0.155 


0.155 


0.155 


0.155 


fuH 


0.051 


0.051 


0.051 


0.051 


0.051 


Power 

Watts 


14 


12.5 


11 


9.5 


8 


Orientation 


Horizontal 


Horizontal 


Horizontal 


Horizontal 


Horizontal 


B9_14 

BN 

3 

0.164 

14 

Horizontal 

B9_125 

BN 

3 

0.164 

12.5 

Horizontal 

B9_11 

BN 

3 

0.164 

11 

Horizontal 

B9_9.5 

BN 

3 

0.164 

9.5 

Horizontal 

B9_8 

BN 

3 

0.164 

8 

Horizontal 

14 

Horizontal 

12.5  ! 

Horizontal 

11 

Horizontal 

9.5 

Horizontal 

8 

Horizontal 

14 

Horizontal 

12.5 

Horizontal 

11 

Horizontal 

9.5 

Horizontal 

8 

Horizontal 

14 

Horizontal 

12.5 

Horizontal 

11 

Horizontal 

9.5 

Horizontal 

8 

Horizontal 

4.8 

Horizontal 

14 

Horizontal 

12.5 

Horizontal 

11 

Horizontal 

9.5 

Horizontal 

8 

Horizontal 
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Run 

Name 

Ceramic 

msmm 

BN 

E2_125 

BN 

E2_11 

BN 

E2_9.5 

BN 

[^■■1 

BN 

E3_14 

BN 

E3_125 

BN 

E3_11 

BN 

E3_9.5 

BN 

E3_8 

BN 

E4_14 

BN 

E4_125 

BN 

E4_11 

BN 

E4_9.5 

BN 

E4_8 

BN 

F1_14 

AIN 

F1_125 

AIN 

F1_11 

AIN 

F1_9.5 

AIN 

Size 

Microns 


Load 

Grams 


0.100 


0.100 


0.100 


0.100 


0.100 


0.150 


0.150 


0.150 


0.150 


0.150 


0.248 


0.248 


0.248 


0.248 


0.248 


0.049 


0.049 


0.049 


0.049 


0.049 


Power 

Watts 


14 


12.5 


11 


9.5 


8 


14 


12.5 


11 


9.5 


8 


14 


12.5 


11 


9.5 


8 


14 


12.5 


11 


9.5 


8 


Orientation 


Horizontal 


Horizontal 


Horizontal 


Horizontal 


Horizontal 


Horizontal 


Horizontal 


Horizontal 


Horizontal 


Horizontal 


Horizontal 


Horizo-^tal 


Horizontal 


Horizontal 


Horizontal 


Horizontal 


Horizontal 


Horizontal 


Horizontal 


Horizontal 


F2_14 

AIN 

44 

0.103 

14 

Horizontal 

F2_125 

AIN 

44 

0.103 

12.5 

Horizontal 

F2_11 

AIN 

44 

0.103 

11 

Horizontal 

F2_9.5 

AIN 

44 

0.103 

9.5 

Horizontal 

F2_8 

AIN 

44 

0.103 

8 

Horizontal 

F3_14 

AIN 

44 

0.207 

14 

Horizontal 

F3_125 

AIN 

44 

0.207 

12.5 

Horizontal 

F3_11 

AIN 

44 

0.207 

11 

Horizontal 

F3_9.5 

AIN 

44 

0.207 

9.5 

Horizontal 

F3_8 

AIN 

44 

0.207 

8 

Horizontal 

Run 

Name 

Ceramic 

Size 

Microns 

Load 

Grams 

Power 

Watts 

Orientation 

F4_14 

AIN 

44 

1.11 

14 

Horizontal 

F4_125 

AIN 

44 

1.11 

12.5 

Horizontal 

F4_11 

AIN 

44 

1.11 

11 

Horizontal 

F4_9.5 

AIN 

44 

1.11 

9.5 

Horizontal 

F4_8 

AIN 

44 

1.11 

8 

Horizontal 

V0_14 

None 

NA 

NA 

14 

Vertical 

V0_125 

None 

NA 

NA 

12.5 

Vertical 

ilSDHii 

None 

NA 

NA 

11 

Vertical 

V0_9.5 

None 

NA 

NA 

9.5 

Vertical 

V0_8 

None 

NA 

NA 

8 

Vertical 

AIN 

44 

5.07 

14 

Vertical 

V3A_125 

AIN 

44 

5.07 

12.5 

Vertical 

SSJBHi 

AIN 

44 

5.07 

11 

Vertical 

V3A_9.5 

AIN 

44 

5.07 

9.5 

Vertical 

V3A_8 

AIN 

44 

5.07 

8 

Vertical 

V4A_14 

AIN 

44 

10.33 

14 

Vertical 

AIN 

44 

10.33 

12.5 

Vertical 

Esm 

AIN 

44 

10.33 

11 

Vertical 

AIN 

44 

10.33 

9.5 

Vertical 

ms^ 

AIN 

44 

10.33 

8 

Vertical 

AIN 

44 

full 

8 

Vertical 

SI2QBI 

BN 

3 

1.95 

14 

Vertical 

mtjm 

BN 

3 

1.95 

12.5 

Vertical 

msoM 

BN 

3 

1.95 

11 

Vertical 

BN 

3 

1.95 

9.5 

Vertical 

iZEOH 

BN 

3 

1.95 

8 

Vertical 

mbom 

BN 

3 

.166 

14 

Vertical 

BN 

3 

.166 

12.5 

Vertical 

BN 

3 

.166 

11 

Vertical 

BN 

3 

.166 

9.5 

Vertical 

SESH 

BN 

3 

.166 

8 

Vertical 

59 


Run 

Name 

Ceramic 

Size 

Microns 

Load 

Grams 

Power 

Watts 

Orientation 

E1B_14 

None 

NA 

NA 

14 

Excitation 

E1B_125 

None 

NA 

NA 

12.5 

Excitation 

E1B_11 

None 

NA 

NA 

11 

Excitation 

E1B_9.5 

None 

NA 

NA 

9.5 

Excitation 

None 

NA 

NA 

8 

Excitation 

E2B_14 

BN 

3 

0.22 

14 

Excitation 

E2B_125 

BN 

3 

0.22 

12.5 

Excitation 

E2B_11 

BN 

3 

0.22 

11 

Excitation 

E2B_9.5 

BN 

3 

0.22 

9.5 

Excitation 

BN 

3 

0.22 

8 

Excitation 

E3B_14 

BN 

3 

1.22 

14 

Excitation 

E3B_125 

BN 

3 

1.22 

12.5 

Excitation 

E3B_11 

BN 

3 

1.22 

11 

Excitation 

E3B_9.5 

BN 

3 

1.22 

9.5 

Excitation 

BN 

3 

1.22 

8 

Excitation 

FC1B_20 

None 

NA 

NA 

20 

Forced 

FC1B_14 

None 

NA 

NA 

14 

Forced 

None 

NA 

NA 

8 

Forced 

FC2B_20 

BN 

3 

1.0 

20 

Forced 

FC2B_14 

BN 

3 

1.0 

14 

Forced 

BN 

3 

1.0 

8 

Forced 

FC3B_20 

BN 

3 

3.66 

20 

Forced 

FC3B_14 

BN 

3 

3.66 

14 

Forced 

FC3B_8 

BN 

3 

3.66 

8 

Forced 
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APPENDIX  C 
SAMPLE  UNCERTAINTY 
CALCULATIONS 
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The  accuracy  of  the  data  in  this  study  was  determined  by  performing  an  uncertainty 
analysis.  For  a  function  comprising  a  number  of  independent  measurements,  F=F(X;,^,  Xg, 
Xc),  the  uncertainty  of  F  was  calculated  as; 

5F  =  [(^5Xa)^  +  CsXb)^  -k  (^5Xc)2] 

If  the  function  F=  K*X*;v*^b*^c»  the  uncertainty  was  determined  by; 

The  properties  of  concern; 


Tp,^=23.9“C 

5Tp,^=  0.25T 

T^=69.5X 

5T^=0.25"C 

Q  *  14.0  watts 

5Q  =  0.25  watts 

p  =  1654  kg/m’ 

6p  =  4.54  kg/m’ 

p=  0.00148  K-' 

5P=  3.46E-6K-' 

Cp= 1 004  joule/kgK 

6Cp=2.4  joule/kgK 

v=4.603E-7mVs 

6v=l  .64E-8m’/s 

a=3.585E-8mVs 

5a=1.26E-10m’/s 

k=0.058W/m*K 

6k=4.30E-5  W/m*K 

L=O.009  m 

5L=  0.00009  m 

A=  0.001 137  m’ 
g*9.807  m/s’ 

6A=1.0E-4m’ 

Uncertainty  calculations  are  for  the  3  micron  BN  light  load  at  maximum  power,  that  is 
mean  power  input  of  14  watts. 
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AT  —  T  s  ~  T  plate 

XT  Q*L 

~  k*A»(AT) 

Nu=  33.89 

With  the  properties  listed  above  8Nu  =2.36  or  6.97  %.  The  Rayleigh  Number  (Ra)  is 
dependent  on  P,  c^,  k,  v,  p  and  a,  all  of  which  vary  with  temperature  for  FC-75  [Ref  21] 

«  _  0.00246 

P  ~  I  .g25-0.00246*T 


Cp  =  (0.2411 +3.7037E- 4  *T) 


k  =  0.065-7.895E-5*T 


u  =  [1 .4074  -  1 .96E  -  2  •  T  +  3.8018E  -  4  *  -  2.731E  -6*T^+ 

8.168E-9*T'‘]*10-^ 

p  =  (1.825-0.00246*1)  *  100 


a  = 


P*Cp 


Thus 


^  =  [(f)'  +  (4f )'  +  (f )'  +  (f )'  +  (f )'  +  (§)'  +  (v)'l 


1/2 


where  Ra  =  ;  Ra  =  9.93E8  and  5Ra  =  8.61E7  or  8.67%. 

kavA 
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appendix  d 

NUMERICAL  ANALYSIS 


PROGRAM 


cxx:cxxcccxm:ccccccccccccxxx:ccccccccccccccccccccccccccccccccccc 

c 

C  AGEhERALPlJRPOSEFCRTTWNPROaiAMFORSOLVINKjTHREE-  $ 

C  DIMENSHWAL  HEAT  TRANSFER  AND  aUID  FLOW  IN  RECTANGULAR  S 
C  COORDINATES  $ 

C 

C  VERSION  CK^PROOIAM  WHICH  CAN  USE  LARCE  ARRAYS  Ftm 
C  FINER  GRIDS  AND  DO  USTEADY  CALCULATIONS 
C  CALCULATIONS  FOR  LARGE  ARRAY  VERSION  ARE  THE  SAME  AS  THOSE 
C  FC»  THE  LOWER  DIMENSION  PROCaL^MEXCPET  THAT 
C  COEFFICIENTS  FOR  F.D.  EQUATIONS  ARE  RECLALCULATED 
C  RATHER  THAN  STORED  AND  REUSED  LATER 

C  THIS  REQUIRES  APPROXIMATELY  25  PERCENT  MORE  TIME  FOR  EXECUTION 
C 

cC  DEVELOPED  BY: 

C  S.B.SATHE 

C  COTE  69ST.  DEPARTMENT  OF  MECHANICAL  ENGINEERING 
C  US.  NAVAL  POSTCKADUATE  SCHOOL 
C  MONTEREY,  CA  93943 
C  TEL:(408)-646-2417 
C  BiTNET  ADDRESS:  5 140P(SNAVPGS 
C 

C  MOTIFIEDFC»  UNSTEADY  OPERATION  BY: 

C  D.  E.  WROBLEWSia 

C  COTE  69WRO,  DEPARTMENT  OF  MECHANICAL  ENGINEERING 
C  U  S.  NAVAL  POSTGRADUATE  SCHOOL 
C  MONTEREY,  CA  93943 
C  TEL:  (408>646-2463 
Relaxabon=0.4 
cDT=l. 

CTnoug|t=0.75 

CEPST=.03 

cRA=I.15e3 

C  uses  logic  fix' oonvogBioe 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

c  LOaCALLSTC*' 
c  INTE<ER*4NOW(14) 
c  CX)MM(»WTm/LSTt*>4CALLJSTC»» 

IXXHCAL  LSOLVELPRINrLKXLSTOT 

COMMON  F(35,25,33,5XP(33,25,35XRH0(35,23,33XGAM(35,25,35), 

1  (XWP5,23,35),AiaX35,25,35),APCM(35,25,35),AP(35,2535), 

2  AIP(35,2545),AIM05  J5,35),AJP(35,25,35)4UM(35,2535) 

COMMON  deIli(35,25,33).deIll0(35a535),e|JH(35,2535X 

3  X(35)XU(35)4(DIF05)4CCV05)4CCVS05Xlk.cjwiIalent, 

4  Y(35XYV(35).YDIF(35XYCV(35XYCVS(35X«nielt.Iprev(35,25,35), 

5  Z(35),ZW(35)2DIF(35),ZCV(35)7CVS(35),ap0(3545,35), 
6YCVR(35),YCVRS(35),ARX(35),ARXI(35)4«UaP(35),apl(35,25,35), 

7  R(35XRMN(35).SX(35XSXMN{35)4CCVI(35)4CCVIP(35), 

8  YCVK35),YCVJP(35XZCVK(35)^CVPa>(35),st(35a5,35) 
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CX»»40NDU(33ja;}5)J)V(35J3,33XDW(33J3.35)JV(35)JVP(35). 

lFX03)JXM(33)JY(33>j^33)JPT(35X<3Tr(35).TOLJD(33.25.35). 

2FZ(35)J=2M(33XVHAT(33^.33).WHAT(35^.35XUOLD(33J3.33) 

CX3MMON/IhDX/RElJO((  13)iJ1^  i  l)J'niMES(  lOX 

lLSOLVE(10XTIMEjyr^Yl^SJWOCON7ERO.TLAST. 

2NF>IFMAX^J«HO>Kl\^aliZUJ^I^I2J^ 

3ISTJSTJCST.rreRXAST. 

4IPREFJPREFJKPREFJ40DE 
COMMON«EADI>WTILE 
CHARACreR*  10  TnLE(  13) 

CXDMMCH4A:^^IlyLSTOPJCAL^^ 

a3MMC»^m4Vl/EPSUJPSV£PSW£PSTJ(X)NV.ITCRl.T0(35.25.35)JENBAU 
1  UO(33a3.33XVO(33a3^3XWO(33a3,35XrrERL 

(XMVfC»4/SOROSMAXSSUM 

ooinniofi/kfoe6(fi)foe(33^,33),yfoioe(33^.3SXzfoice(33^S.3S) 

(XWNO4«D0EF/F1j0WJDIFFA00F 

DlMENSK»4U(33a3,35).V(33a3.35XW(35J13.35)J>C(33J5.35) 

DIMENSK»<T(33^3.33) 

EQUIVALENCE(F(1.1.1,1XU(1.1.1)).(F(1.1.U).V(1.1.1)). 

1  (F(1.1,1.3XW(1.1,1)X(F(1.1.1,4XPC(1.1.1)) 

2  .(F(1.1.1.3XT(1,1.1)) 

afien(16jGle>‘'pooanf) 

q)en(20Jle^iftiogO 

ISTOIM) 

CALL  MESH 
CALLGEOMET 
CALL  BEGIN 
lOCALLVARRHO 
19999  PORMATC  LSTOP-  ’42) 

CALLBNDRY 
CALLPRTOUT 
CALLINCRV 
IF(LSTOP)'TTffiN 
WRrrE(M9999)ISTOP 
STOP 
ENIHF 
CALLCOEFF 
GOTO  10 
END 

(XCCCXXCXrCCCCCCCCCCCCCCCCCXXXXXXXCCCCCCCCXXXXrCCCCCCCCCCC 
SUBROUTINE  FROnL 

ccccccxxxxccxxxxxxxxcccccccccccxxxxxxmmxcccccccccccccc 

CX3MM»^CX3EF/nX)WJHFFACOF 

ACOF-DBFF 
quint* ACOFDIFF 
c  flow^. 

1F(FLOW£Q.O.)  RETURN 
TEMPKMFF-ABS(FLOW)*0. 1 
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ACX3F-0. 

IFOEMPLEO.)  RETURN 

TEKff^TEMP/DIFF 

ACOF-DIFPTEMP**5 

RETURN 

END 

(xmxr(xxxcx:(xccxxxxx:ccxxcxmxxxxxccccccccccx:cccx:cxxx:ccc 

SUBROUTINE  H^lA 

(xccccccrc(xccccccccccc<xccccc<xccaxxxrccx:ccccccccccccccc 

LOGICAL  LSOLVEJJRINTimiCLSTOP 

COMMON  F(35^.35,5XP(35^,35)JIHO(35^  J5XGAM(35^.35). 

1  CON(35a5,35),AKP(35^^5)u\KM(35^^5XAP(35^.35). 

2  AIP(35J5,35XAIM(35a5;35XAJP(35^,35XAJM(35^.35) 

COMMON  ddh(35^.35).deUiO(35^;}5),epa(35^.35), 

3  X(35)^35);CDIF05)^CV(35);CCVS(35Xllwi>,afaieBl, 
4Y(35XYV(35XYDIF{35).YCV{35XYCVS(35).tineft.lpitv(35^.35X 

5  ZP5XZW(35XZDlF(35XZCV(35)7CVS(35Xap0(35;i5,35). 

6  YCVR(35).YCVRS(35).ARX(35XARXJ(35).ARXIP(35)^1(35^.35). 

7  R(33XRMN(35XSX(35XSXMN(35XXCVK35XXCV1P(35). 

8  YCVJ(35XYCVIP(35).ZC\^35),ZCVKP(35Xs«35^.35) 

G(mfC»4  DU(35;55.35XDV(35^.35XDW05;i5^5XFV(35)JFVP(35). 

1  FX(35)JXM(35)JY(35)JYM(35XPT(35).Qr(35XTOLD(35^.35). 

2  EZpS)JFZM(33XVHAT(35^,35XWHAT(35^^5XUaLD(35^,35) 
COMMON/INDXmAX(13)iJWr(13).LBLK(l  1XNI1MES(10), 
lLSOLVE(10XTIMEJ)T^YL^SJ«Oa»i.ZERO.TlAST. 
2NFJ4FMAX^>IRHO>IGAM^lJJiJJMl>122^J41J422^^ 
3ISTJSTJCSTJIERJLAST. 

4IFREFJFREFja>REFAI0raE 

DIMENS04  D(33XVAR(3SXVARM(35XVARFpSXPHlBAR(35) 
OQfniiioii/fbice/'xfbioe(3S^^SXyt)ice(3S^,33Xzfxo^ 

COMMONMEADINnrnLE 
CHARACTER*!©  TrTLE(13) 

C********************************************************* 

CCCCCCCXXXXXXXXXXXCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

OOMMON«DOEF/FLOWJ)IFF>VCOF 

ISTFHST-l 

JSTF«JST.l 

KSTF=KSr-l 

m*u+isr 

rrz-u+BT 

JTl^C+JST 

JT7-M3+JST 

KT1-N2+KST 

KTW434KST 

C - 

NTSSS-NTTMESCNF) 

qxmt*^,NrSSS 

D09991>n^ljnSSS 

D0391N>^J4F 
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c 

IF(.NOrr.LBLK(NF))CJO  TO  60 
60  CONTINUE 

COMMENCE  TOMA  LINE  BY  UNE  SWEEPS  FOR  SOLUTION 
DO90K-KSTW 
DO90WSTM2 

Fr(isn><). 

QT(ISTT)-F(ISTFJJC.N) 

DO70I-ISTJ2 

CCNOM-AP(UXHTa-l)*AIM(UJC) 

Pr(D-AIP(UJCVDENOM 

TEMFKX)N(UJCKA;P(UJC)*F(I>lJCN><-AMUX)*F(U-liLN) 

1  +AKP(UX)*F(UX+I.N)*-AKM(UX)*F(UX-1X) 
QT(IKIEMF^AIM(UXWa-l)yDE^ 

70CC»mNUE 

DO80n=ISTX2 

i-m-n 

80  F(UXNH^a+lJJCN)*PT(I>+<3T(I) 

90cca>mNUE 

c - 

DOIWKK’^CSTXB 
K>Kn-KK 
DO  190  JJ-JSTJM3 
J-m-JJ 
PT(isn>o. 

Qfr(ISTI>F(ISTFJJtN) 

DO  170I-ISTX2 

DENOM-AP(UX>PT(I-irAIM(UX) 

Pr(D-AIP(LJX)/DENOM 

TEMPKXDN(UX>^AJP(UXfF(I>lXJ<>^-AJM(UXrF(U-IXN) 

1  +AKP(UJO*F(UX+lXHAKM(UXrF(UX-IX) 
Qr(IKIEMPfAIMaJX)*Qra-l)yDENOM 
ITOCWITNUE 
DOiaoiHsrjLz 
MTi-n 

180  F(UJLNH'a+UXJ'im(IH3T(I) 

190OOFniNUE 

C - 

DO290HSTX2 

DO290K-KSrX2 

PT(JSIE)H). 

QTCJSTFHWSTFJLN) 

DO270>-JSrXQ 

DENQM-AF(UX>fTCJ>l)*AJM(UX) 

FT(J>-A)P(UXyDENOM 

TEMIHDON(UX>fAKFXUX)*F(UX+IX)^AKM(UX)*F(UX*W 

1  +AiP(ujo^(^uxx)^-AiMajx)*T^a-ux^ 

QT(JKIEMP4>A1M(UJ0*QT(J-1)VDEN0M 

TJOCOtmSUE 

DO280JI-^JM2 
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i-jri-jj 

280  F(UJCJ4Hf(l>l  ja^PT(J>+<7T(J) 

290CX»mNUE 

C' - - - 

DO390Ih‘ISTiJ 

I-IT2-n 

DO390KK=»KST.N3 

K=4CT2-KK 

Pr(JSTPH). 

QT(JSTTHf(USTTJCN) 

DO370>=JSTJ^ 

DEhK)M“AP(UJC)^(J-l)*AJM(UJC) 

PT(J)-AJP(UJ^yDENOM 

TEMP<X»J(UX)+AKP(UJC)*F(UJC+l>I)+-AKM(UJC)*FaJJC-l>l) 
I  +AowjcrFa+uja<>+AiM(UJC)*Fa-ujcN) 
Qrr(J)-<TCM^-AJM(UX)*QT(J-l)yDEh^ 

370CX»rnNUE 

DO380JJ=JSTX12 

J-JTl-JJ 

380  F(UXN>-F(U+IXJ'0*PT(J>+<3T(J) 

390CCWTINLE 

- — - 

DO490WSTXI2 

DO490I»ISTiJ 

PT(KSTFH). 

Qrr(KSTF)-F(UXSTFJ4) 

D0470K-KSr>[2 

DEN0M-AP(UXHT(K-1)*AKM(UX) 

Pr(K)-AKP(UXyDENOM 

TCMFKX»4(UJK)+AIP(UX)*Fa+lJXJJ)<-AlM(UJO^ 

I  +AJP(UJC)*F(W^lXW-AJMaJJO*F(U’lXJ^ 

QT(KKreMPfAKM(UX)*Qr(K-l)>I>ENOM 

470CC»4nNUE 

DO480KK»KST^ 

K-icri4ac 

480  F(UXJ4)-F(UX+i  WPT(K>K3T(K) 

490a»mNUE 

c- - - - 

DO590JJ-JSTXO 

F-JR-XI 

DOSWIt-ISTiJ 

i-fn-n 

PT(KSIFH). 

Qrr(ICSIF)^(UXSTF;4) 

DO570K-*OTJC 

DENO!^AP(UJO^>TXK-irAia^X) 
FT(K>-AKP(UXyDENOM 
TCMP<X)hWX>+^AIP(UX)*F(H-UXJ<)fAlM(U 
I  +AJP(UX)T^(U+iXJ<)+-AJM(UJC)*Fa>lXJ^ 
QTaCHTEMP^AKM(yjCr<^ 

STOCXWITNUE 


69 


D05»KKH(STJ42 

K-KTl-KK 

580  F(UJCN)-f(UJC+l  J4)*PT(K>K3T{K) 

590  CONTINUE 

C - 

391  CONTINUE 

999  CONTINUE 
call  reset 
return 
end 
C 

Q**************..**********,***.,..,******..**.* 

subroutine  RESET 

e*********************************************** 

LOGICAL  LSOLVEJLPRINTJLBLKXSTOP 

COMMCWF(35a5.35,5)J»(35J5.35)JRHO(35a5.35).GAM(35.25.35). 

lCON(35J5,35),AKP(35J5.35).AKM(35J5^5)uM>(35^.35). 

2  AIP(35^5.35)jaM(35J»5,35)>UP(35^.35).AJM(35^.35) 

COMMCX^I  deili(35^.35).(UliO(35^,35).q)Si(35aS.35). 

3  X(35XW(35)^IF(35XXCV(35)^CVS(33Xtk.cp.alaient, 
4Y(35).YV(35).YDIF(35).YCV(35).YCVS(35).tinelt,l|)iev(35a5,35X 

5  Z(35UW(35XZDIFa5),ZCV(35)7CVS(35Xa|j0(35^.35). 

6  YCVR(35).YC\^5).ARX<35).ARXJ(35XARXJP(35XapU35^.35). 

7  R(35)JtMN(35XSX05XSXMN(35);CCVI(35XXCVIP(35X 

8  YCVK35XYCVJP05)JX:VK(35)jZnia\3^ 

COMMC»I  DU(35^  J5)J)V(35a5,35XDW(35^  J5)JV(35XFVP(35). 

1  FX(35)J^35)JYO5)fYM05)jTip5XQTO5XTOIJX35^,35X 

2  FZ(35XFZM(35XVHAT(35^;J5),WHAT(35^35XUOLD(35i535) 
CXX^1M(:»l/INDX/RELAX(13).LPRINni3)JA^ 
lLSOLV^10).T1MEXrT^YL,ZL.SJtHOCC8i2m).lL^ 
2NFJlFMAXJ4PJ9RHOJ«jAMJ-lJ2^3>Il>OJ^Jll,f^ 
3ISTJSTJCSTJIERMST, 

4IPREFJPREFJKFREFMX3E 

COMMONMEADINariLE 

CHARACIER*10  HILBCia) 

COMMONOmASTOPJCALLJSlOP 

COMMOWOONVl/EPSUfPSVJEPSW,EPST4CONVJlERl.TO(35^35XENBAU 
1  UO(35^45).VO(35^35).WO(35^.35XnERL 

COMKON«ESIDmiiAX(13XIRESID 
oooiinai«KW^cfiitaP5^45)jfac^ 

COMMON/SQROSMAXSSUM 

(XIMMONCCEF/FLOWJMIT^U^ 

DIMENSION  U(35^33).V(35^35).W(35;2545XPC(35^.35) 

DIMENSION  T(35^^5) 

EQUIVALENCE(F(l,l,I.l),U(l,ia)X(F(LLU).V(LI.l)), 

1  (F(Ll.l,3XW(l.l,l)).(F(Ll,l,4XPqLLl)) 

2  .(F(1.1.1AT(1,1.1)) 

D0400K=-2^ 

D0400>>2>Q 
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DO  400 1-2X2 
CXDN(UJC)-0. 

AP(UJ0=O 

400CXKnNl)E 

RETURN 

END 

C 

crccccccccccccccccccccccccccccccxxrccccccccccccccccccccccc 

SUBROUTINE  FORM 

ccccccc<xccx:ccccccc<xcc(x:cccxxxxccx:ccccccccccccccccccccccc 

LOGICAL  LSOLVEXPRINTXBLKXSTOP 
COMMC»4F(35.25.35,5)J^35.25,35)JIHO(35.25.35).GAM(35.25.35). 
lCC»i(35^5.35).AKP(35^.35),AKM(35^.35)^35.25.35). 
2AIP(35,25.35XAIM(35,25.35)AIP(35X3.35)w^JM(35.25.35) 

COMMON  <fclh(35^3.35),cldh0(35a5.35),qpsi(35.25.35), 

3  X(35)^(35)^IF(35)J(CV(35)^C\^35Xflc.cp.alaient. 
4Y(35XYV(35).YDIF(35).YCV(35).YCVS(35).troAtprev(35^.35). 

5  Z(35)^(35)^lF(35).ZCV(35)7CVS(35),ap0(35^.35). 

6  YCVR(35),YCVRS(35).ARX(35).ARXJ(35).ARXJP(35)jpl(35^.35). 

7  R(35XRMN(35XSX(35XSXMN(35);CCVK35XXCV1P(35), 

8  YCYK35).YCVJP(35)2CVK(35)7CVKP(35)^35^,35) 

COMMCW  DU05X5.35)X>V(35a5,35U)W(35;i535)J^(35)XVP(35X 

1  FX(35)JTCM(35)JY(35)J^35)JT(35).<F(35XTC)LD(35^.35), 

2  FZ(35XFZM(35).VHAT(35^.35),WHAT(35^35),U0LD(35^.35) 
C(XvIMCX4^X/RELAX(13)XPRINT(13)XBL^  1).NTIMES(10X 
lLSOLVE(10XTIMEXn’;CL.YL7LARHOCON.2ERO.TLAST. 
2Mf>IFMAXJIPJ«HOJ«AMXli2X3Ml>I2J^J41>CXD. 
3ISTJSTJCSTjrERXAST, 

4IPREFJFREFja>REFMODE 

C(Xk4MC»I/HEADIN/TTILE 

CHARACTER*  10  TnLE(l3) 
a»IMa4/CmiAST»4CAlXJST^ 

C0MM0N«D0NV1/EPSUEPSVEPSWEPST4C0NVJIER1,T[)(35^35XENBAX 

1  UO(35XS^5),V[)(35X5^5).WO(35^.35>4IERL 

COMMC»4/RESIEVRMAX(13)4RESID 

COMMC»«SOROSMAXSSUM 

OOMMONCOEF/FLOWZMFFACOF 

DIMENSION  U(33^35),V(35^;35),W(35^^5XPC(35^.35) 

ooinraoDlbR»^xfxcc05^3S)jif>iOG0945^S),z&ice^ 

DIMENSION  T(33^35) 

EQUIVALENCE(F(l,l.l.l),U(l,l.l)X(F(l,l.UXV(l.l.l)). 

1  (F(l.l.UXW(l,l,l)X(F(l,l.l,4XPC(l.l.l)) 

2  ,(F(1,I,1AT(1,1.I)) 

1  F(3RMAT(13X,XXIMFl7rA'IK)N  IN  CARTESIAN  COORDINATES') 

2  FC»MAT(lSX,t]OMFUTATKX4  FOl  AXISYMMEIRIC  STTUATICK) 
3FORMAT(13X,'COMPin'ATION  IN  POLAR  COORDINATES') 

4  PC»MAT(14X40(1H*X//) 
itHBn 

end 

C 

COME  HERE  TO  CALCULATE  CiRII)S  SPECIHCATION 
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o  o 


subroubne  GEON/ET 


LOGICAL  LSOLVEiJRim'XBLICXSTOP 
a)MMCWF(35.25,35.5)J»(35J5.35)JRHO(35.25.35).GAM(35.25.35). 
lCON(35,25.35).AKP(35J5.35).AKM(35.25.35)^35.25.35). 
2AIP(35.25.35).AIM(35.23.33).AJP(33J3.35)AIM(35^3.35) 
(XM^Nddh(33,23.33).deUiO(33J3,33).epsj(35^3.33). 

3  X(33)jaj(33)^IF(33UCV(33)^CVS(33),tk.cpjilatenL 
4Y(33).YV(33).YDIF(33),YCV(33),YCVS(33).lindLlpitv(33^3.33). 
3Z(33)jZW(33)7DIF(33)2CV(33).ZCVS(33)^33^.33). 
6YCVR(33).YC:VRS(33),ARX(33).ARXJ(33XARXJP(33)^1(33.23.33). 

7  IU33)JIMN(33).SX(33XSXMN(33)^CVI(33)J<CVIP(33). 

8  YCVJ(33),YCVJP(33)  jSCVK(33).ZCVKP(33)^33.23.33) 
CX)MMONDU(33a5.33)i)V(33.23.33)JDW(33a5.33XFV(33)JVP(35), 
lFX(33)J^<^^33)JY(33)J^lM33)J^’(33X<F(33),TOLD(33^.33). 
2FZ(33)J7H33),VHAT(33^,33).WHAT(33J5J5).L)0LD(33J3.33) 
(XM^I(M^/INDX/RELAX(  13UJ>W^f^  13)J-BLK(  1  l)a^nMES(  10). 
lLSOLVE(10),TIMEJDTJ<L,YLiLSJlHO(XW2ERD.TLACT 
2NF>IFMAX>lP.NW«J4GAMI.lJL2iJMlMLM3>IIJ^ 
3ISTJSTJCSTJTERXAST, 

4IPREFJPREFja>REF>«X)E 

CX3MMON«EADINrrrrL£ 

CHARACTCR*l0TraE(l3) 

OIMMONOmASTOPJCALLJST^ 

OO^«^iONmWl/EPSU£PSVfPSW£PST4CX)NVJTCRi.T0(35a5,35XENBALs 

1  U0(33J3.33).VD(33.23.33XWO(35a5.35XrreRL 

CXX^fMC»4/RESID/RMAX(13)JRES^ 

ooinroon/lbrceWbrce03.23,33),ylbrce(33.23,33).2fi)rce(33.23.33) 

COMMON/SORC/SMAXSSUM 

OCM^4C»4/COEF/FLOW4>IFF.ACOF 

DIMENSION  0(35.2335).  V(33a3,33),W(35.25J5XPC(35.2535) 

DIMENSION  T(333535) 

EQUIVALENCE(F(I.I.I,IXU(I.LI)X(F(I.I.UXV(1.I.1)X 

1  (F(I.I,I.3XW(I.1,I)X(F(I.1,I.4XPC(I.1.1)) 

2  .(F(1.I.I,3XT(1.I.I)) 

1  FOW^T(l5XXDC»«VrATICW  IN  CAKIESIAN  COORDINATES') 

2  FORMAT(I5X.'COMIVrATION  FOR  AXISYMMEITUC  SnUATTON^ 
3FORMAT(I3X.'COMPlJrATION  IN  POLAR  COORDINATES^ 

4  FORMAT(14X,40(Un//) 

L2-LM 

L3^-I 

M2>M1-1 

M3-M2-1 

N2-N1-1 

N3«N2-I 

X(I)»XU(2) 

D05I=2X2 

5  X(I)-0.3*(Xlia+I>^-XU(D) 

Xi(Ll)^XU(LI) 
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Y(1>-YV(2) 

DO10J»2M2 

10  Y(J)=0.5*(YV(J+  l)+YV(J)) 
Y(M1)=YV(M1) 

Z(1)=ZW(2) 

D07K»2J^2 

7  Z(K>=0.5*(ZW(K+IKZW(K)) 
Z(N1)=ZW(N1) 

C - 

DO  15 1=2X1 
15XDIF(D=X(D-Xa-1) 

DO  18 1=2X2 
18  XCV(D=XUa+l>-XU(D 
EX)  20 1=3X2 
20XCVS(I>=XDIF(I) 
XCVS(3)=XCVS(3>4-XDIF(2) 
XCVS(U>=XCVS(UHXDIF(L1) 
DO  22 1=3X3 
XCVI(I)=0.5*XCV(I) 

22  XCVIP(I)=XCV1(D 
XCVIP(2)=XCV(2) 
XCVI(U)=XCV(L2) 
D0175K=2J^1 
175  ZDIF(K)=Z(K)-Z(K-1) 
D0178K=ZN2 
178  ZCV(K)f=ZW(K+1VZW(K) 
DO270K=3.N2 
2702CVS(K)=ZDIFaQ 
ZCVS(3)»ZCVS(3>+^ZDIF(2) 
ZCVS(N2)=ZCVS(N2)+^ZDIF(N1) 
D0272K=3J43 
ZCVK(K)=0.5*ZCV(K) 

272  2CVKPaC)»ZCVK(K) 
ZCVKP(2>=ZCV(2) 
ZCVK(N2)=ZCV(N2) 
D035J=2M1 
35YD1F(J)=Y(J)-Y(J-1) 
D040J»2>f2 
40  YCV(J)-YV()<-1)-YV(J) 
D045J-3M2 
43YCVS(J)-YD1F(J) 
YCVS(3)-YCVS(3)+YDIF(2) 
YCVS(M2)=YCVS(NC)fYDIF(Ml) 
D0277J-3AO 
YCVJ(J)-0.5*YCV(J) 
277YCV3P(J)-YCV3(J) 
YCVJP(2)=YCVa) 
YCVJ(M2)=YCV(M2) 
IF(MODEJ4E.1)GOT055 
D052MM1 
RMN(J)=1.0 


52  R(J)=l.O 
GO  TO  56 
55DO50J=2Ml 
50R(J)=R(J-1>+-YDIF(J) 

RMN(2>=R(1) 

DO60J^3M2 

60  RMN(J)=RMN(J-I>+-YCV(M) 
RMN(M1)=R(M1) 

56  CONTINUE 
D057J=1M1 
SX(J)=1. 

SXMN(J)=1. 

IF(MODE.NE.3)  GO  TO  5^? 

SX(J)=R(J) 

IF(J.NE.  1)  SXMN(J>=RMN(J) 

57  CONTINUE 
D062J=2>12 
YCVR(J>=R(J)*YCV(J) 

ARX(J)=YCVR(J) 

IF(M0DE.NE.3)  GO  TO  62 
ARX(J)=YCV(J) 

62  CONTINUE 
D064J=4>D 

6t  YCVRS(J>=0.5*(R(J>fR(J-l))*YDIF(J) 
YCVRS(3)=0,5*(R(3>fR(l))*YCVS(3) 
YCVRS(M2)N).5*(R(MIHR(NC))*YCVS(M2) 
1F(M0DE.NE.2)  GO  TO  67 
D065J=3M3 

ARXJ(J)^.25*(l.+RMN(jyR(J)rAPX(J) 

65  ARXJP(J)=ARX(J)-ARXJ(J) 

GO  TO  68 
67D066J=3M3 
ARXI(J)=0.5*ARX(J) 

66ARXJP(J)=ARXJ(J) 

68ARX)P(2>»ARX(2) 

ARX}(M2)=ARX(M2) 

DO70J=3M3 

FV(J>=ARXJP(jyARX(J) 

70FVP(J)»l.-FV'(J) 

D085^=3J2 

FX(D>0.3*XCVa*iyXDIF(D 

85FXM(I>*l.-FX(D 

FX(2)=0. 

FXM(2)-1. 

FX(L1)=1. 

FXMCLl)^. 

DO90^3>12 

FY(J)FK).5*YCV(J-iyYDIF(J) 

90FYM(J)-1.J^(J) 

FY(2H>. 

FYM(2)-1. 
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FY(M1>=1. 

FYM(M1H). 

D087K=3,N2 

FZ(K)=0.5*ZCV(K-iyZDIF(K) 

87FZM(K)=1.-FZ(K) 

FZ(2>>0. 

FZM(2>=1. 

FZ(Nl)=l. 

FZM(N1)=0. 

COHAP.U.VJIHOK:  and  P  arrays  are  INmALIZED  HERE 
D095K=l.Nl 
D095J=1JV11 
00951=1^1 
U(UJCH). 

VaJJCH). 

W(UJC)=0. 

DU{UX>=0. 

DVftJJCH). 

DWdJXH). 

CX)N(UJC)=0. 

APOJXH). 

RHO(UJK)=RHOCON 

P(UJC>=0. 

AIP(UXH). 

AJP(UJC)=0. 

AKP(UJC)=0. 

AIM(UJC>=0. 

AJM(UJC)=0. 

AKMOJJCH). 

PC(UJC)=0. 

00N(UJC)=0. 

95CX»mNUE 
IF(MCCe£Q.  1)  PRINT  1 
IF(MODE£Q.2)PR]Nr2 
lF(MCX)EEQ.3)PRINr3 
PRINT4 
RETURN 
end 
C 

CX)ME  HERE  TO  CALCUALTE  OOEITICIENTS  FOR  FINHE  DIFF.  B^S. 

fy*0**00m***********0***m**00****00***^*0*0*00***0*0m***** 

suixouliiieCX^FF 

LXX3CAL  LSCX.VE.LFRim'JLfLI^^ 

CXDMMON  F05^.35.3)TO3^.35)J«0(35^35).GAM(35^,35), 

1  CON(35^;J5)>Kia»(35^,35),AKM(35a545)>tf(3^^^ 

2  AIP(35^,35).AIM(35^.35)^35^.35).AJM(35^;}5) 

CX)MMON  delh(3S^,3S),(ldbO(3S^4S),e|]^S^3SX 

3  X(35XXUP5)^ff(33XXCVa5)^CVS(35Xlk,qwlam 

4  Y(35).YV(35).YDBF(35).YCV(35).YCVS(35).tnidt.tpiev(35^,35), 

5  2(35)^(35)^IF(35XZCV(35)jZCVS(35)^5^,35X 
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6YCVR(35).YCVRS(35XARX(35).ARXJ(35).ARXJP(35).^1(35.25.35). 

7  R(35)JIMN(35),SX(35XSXMN(35)JCCVI(35)J(CV1P(35). 

8  YCVK35).YCVJP(35)vZC\^35)iX^TCP(35).sU35.25^ 

CX)MMC»4DU(35^5,35)J)V(35J5,35)J>W(35^5.35)JV(35)JVP(35). 
lFX(35)J3a4(35)JT(35)J^35)JPT(35).Cn'(35).TOL^ 
2FZ(35)JFZM(35).VHAT(35^.35),WHAT(35i5.35),UOLD(35.25.35) 
COM^«»W^roX/REl^l3)iJ^U^  IXhmMESdO). 

1LS0LVE(  10).TIMEJ)TJCUYl^S  JlHOCX>NimO.T^^ 
2NF.NFMAXNP,NRHO.NGAMX  liJJ-3  Ail>12>D.N  I.N2.N3. 
3ISTJSTJCST,nERMST. 

4IPREFJPREFJCPREFMC»E 

(X)MMON«EADINnrnLE 

CHARACTER*  10  TnLE(l3) 

(X»MN«»40rilASTX»»JCAlJUISTC^ 

<X)MMON/CO^A^l/EPSUEPSVEPSWEPST.ICONV.r^ERl.TO(35.25.35)ENBA^ 

I  U0(35^,35).VD(35.25,35).W0(35.25.35)J1ERL 

C0MN«:»s</RESID«RMAX(I3),IRESro 

ooiiunoti^bioe^xfiMce(3S.25,35),yfoioe(35JS.3S).;;fi>roe(35.2S,3S) 

(XIMMON/SORC/SMAXSSUM 

<XlMMON/CX)EF/njOWX>nT.^ 

DIMENSICWU(35.25.35),V(35.25.35),W(35^5,35)JC(35.25.35) 

DIMENSKW  7(3505,35) 

EQUIVALENCE(F(l,l,l.l),U(l,l.l)).(F(l.l.UXV(l.l,I)). 

1  (F(1.1.1.3),W(l,l.l)X(F(I,l,l,4XPqi.I.l)) 

2  .(F(l.l.l,5).T(l,l,l)) 

C 

C  STARTOPnERATTONUXIPFOREACHTIMESTEP 
C 

«**««««««««««««« 

99  CC»mNUE 
ICALL-1 
C 

C  CALL  UCOF,  WCOF,  AND  WCC*^  TO  FIND  UHAT,  VHAT.WHAT 
CALLUOOF 
CALLVCX)F 
CALL  WOOF 
C 

C  CALLPOOFTOCALCULATEPRESSUREEQN.COEF. 

IF(.NOTiSOLVE(NF))  GO  TO  500 
CALLPCOF 

IF(nERIJE.I)QOTO409 

D0408IC>0,N2 

DO408J-2>Q 

DO406MJL2 

AF(UJQ*AP(UXyRELAX(NP) 

C0^WJCKX)N(UJQH1.-RELAX(NP))*AP(UJ^^ 

4080CX^ITNUE 
409  (CONTINUE 
NF-4 

CALLTDMA 
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c 

C  CALCULATE  U* 

C 

ICAa-2 

CALLUOOF 

NF-1 

IST»3 

JST=2 

KST*2 

D0413K=2J« 

D0413J=2>12 
00413  I=3i^ 

413  C0N(UXKX)N(UJC>4OU(UX)*AP(UJC)*(Pa-UJCVP(UX)) 

C 

C  SAVE  (XJ)  VALUES  OTUFC»CONVER(ENCE  CHECKING 
C 

007041=10.1 

00704J=1M1 

00704K=1.N1 

UOLO(UJC>=U(UJC) 

704  CCM^TINUE 
C 

CALL  TOMA 
C 

C  CALCULATE  V* 

C 

CAaVCOF 

NF=2 

IST=2 

JST=3 

KSr-2 

00513K=2J^2 
005131=2X2 
00513  J=3XC 

CON(UJCKX)N(UOC>H3V(UJCrAP(UJC)*(P(U-lJC)WJC)) 

513oontiniK 

C  SAVE  CU)  VALUES  OF  VFC»  CONVERGENCE  CHECKING  CIEMPSTCRE  IN  VHAT) 
C 

007141=1X1 

007143=1X11 

00714K=lXiI 

VHAT(UK)=V(UK) 

714  CONTINUE 
CALLTI»4A 
C 

C  CALCULATE  W* 

C 

CALLWCC*? 

NF=3 

IST-2 

jsr=2 


77 


KST*3 

D0523>“2JVC 
DO  5231=24-2 
D0523K»3,N2 

(X)hKUJC)=CON(U4CH3W(U4CrAP(U4C)*(P(UX-l>WJ^^^ 

523  oontiiiue 

C 

C  SAVE  OLD  VALUES  OTWFC»CON\^(2NC:ECHEC3aNG  (TEMP  STCKE  IN  WHAT) 
C 

DO  724 1=14-1 
D0724j=lMl 
D0724K=LN1 
WHATaJX)=WaJ4C) 

724  CONTINUE 
CALLTDMA 
C 

COEmCENTS  FOR  THE  PRESSURE  (DORRECTTON  EQUATION - 

C 

CALLPCXF 

NF=4 

E(.Nar4-SCK-VE(NF))  GO  TO  500 
IST=2 

jsr»2 

KST=2 

CAULDOTUS 

SMAX=«. 

SSUM-0. 

C  WRnE(V)'P;(DOEFF2' 

DO410K=2,N2 
DO410>=24V12 
DO  410 1=2X2 

Va-=YCV(J)*XCV(D*ZCV(K) 

410  CCM4(U4C><X»4(U40*V0L 
DO701K=24^2 
DO701J=2>12 

(X>N(2JJCHlHO(U4CrU(2J4C)*YCV(J)*ZCV(K) 

701  CCWITNUE 
DO  702 1=2X2 
DO702K=2Xl2 

<X)N(lXK)=RHO(L  lXrV(L24C)*XCV(I)*ZCV(K) 

702  CXM^nNUE 
DO  153 1=24-2 
DO  1534=2X12 

CON(U;tHlHO(U.l)*W(UJ)*XCV(I)*YCV(J) 

153  CWITNUE 
C 

c  WRnE(**)’P:COEFF3' 

D035IK-2XC 
DO  351 4=24^ 

DO  351 1=2X2 
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AREA-YCV(J)*ZCV(K) 

AW«>AREA*(FXa+irRHOa+UJC>+FXMa-*-l)*RHO(UJC)) 

FLOW^ARHO*Ua+UJC) 

(X>N(UJCHX)N(UJC>FLOW 

(X>Na+U4C)KX)N(I+UJCHljOW 

C 

AREA-XCV(D*ZCV(K) 

ARHC>=AREA*(FY(J+irRHO(U+lJC>fFYM(J+l)*RHO(UX)) 

FL0W=ARH0*V(U+1JC) 

OON(UJC)-CON(UJCH=LOW 

a»4(l>lJCHX)N(I>lJC>+FL0W 

C 

AREA=XCV(I)*YCV(J) 

ARHC>=AREA*(FZ(K+l)*RHO(UJC+l)+^FZM(K+irRHO(UX)) 

FLOW=ARHO*W(UJC+l) 

CX)N(UJOKXDNaJJC>njOW 

(X)N(UJC+lHajOW 

PC(UJCH). 

351  CXDmiNUE 
D0352I=2i^ 

D0352J=2>I2 

D0352K=2,N2 

SMAX>AMAXl(SMAX>^a)N(UX))) 

SSUM»=SSUM+OON(UJC) 

352  C(»mNUE 
CALL  TOMA 

C 

(XIMEHERF  TO  CORRECT  TOE  VELOCmES - 

C 

C  WRnE(**)'VELC(»REC* 

D0511K=2J42 

D0511J=2M2 

D0511I=2JJ 

IF(LNE.2)  U(UJCHJ(U;C)fIXJ(UJC)*(PCa-UJC>«UJJC^^ 
IF(J.NE2)  V(UJC>*V(U4C>+DV(LJ^)*(PC(IJ-1X>«XUJC)) 
IF(K.NE2)  W(U4C)*W(UJC)+DW(UJC)*(PC(UJC.1H>C(UJC)) 
511COmTNUE 
500CCKnNUE 
C  WRnE(**)'TOOEFFl  ’ 

C 

COEmCIENTS  POROTOEREQUATO»1S - 

C 

CALL  RESET 
1ST=2 
JST-2 
KST»2 

D0600NF-5>IFMAX 
IF(.NOTLSOLVE(NF))  GO  TO  600 
CALLDIFFUS 
REL-1.-RELAX(NF) 

DO  601 1-2X2 
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DO601>-2Af2 

DO601K-2J« 

ooEFnoEhrrs  av^-est  and  aeast 
AREA-  YCV(J)*ZCV(K) 

FLOW-AREA*Ua+UJC)*(FXa+l)*RHOa+UJC>H=XMa+l)*RHOaJJK)) 

DlFF-AREA*2.*GAM(UX)*GAM(l+UJCy(XCV(I)*GAMa+lJX)+ 

1  XCVa+irGAM(UX>f  1.0E-2fi) 

CALLPROFIL 

IF(ABS(FLOW)i-T.  lE-20)  FL0W=0. 
AIMa+lJJK)-ACOT+AMAXl(ZEROjTLjOW) 
AIP(UJC)-AIM(l+UjCVFLOW 
COEmCEhnS  ANCWra  AND  ASOUra 
AREA-  XCV(I)*ZCV(K) 

FlX)W-AREA*V(U+lJC)*(FY(J+l)n«0(U+lX>+FyM(J+l)*l^^ 

DIFF=AREA*2.*GAM(UJC)*GAM(I^lXy(YCV(J)*GAM(I>l^ 

I  YCV(J+l)*GAM(UX)+-l.(E-20) 

CALLPROFIL 

AJM(U+lX)=A(Xf+AMAXl(ZEROjajOW) 
AJP(UX>=AJM(U+1JC>FIjOW 
COEFFICIENTS  AOUT  AND  AINTO 
AREA-  YCV(J)*XCV(D 

FLOW-AREA*W(UJC+I)*(FZ(K+l)*RHO(UJC+l>+fZM(K+l)*RHO(UX)) 
DIFF=AREA*2.*GAM(UJC)*GAM(LJJC+iy(2CV(K)*GAM(UJC+l>^- 
1  ZCV(K+l)*GAM(UJC>+-l.(E-20) 

CALLPROFIL 

IF(ABS(FLOW)JLT.  lE-20)  FLOW=<). 
AKM(UJC+l)“ACOF+AMAXl(ZEROJ=LOW) 
AKP(UJC)»AKM(UJC+1>FIjOW 
6oia»mNUE 
C  WR1TE(V)'TC0EFF2  ' 

DO610J-2M2 

DO610K-2J42 

COEFFICIENTS  AWEST  AND  AEAST 
AREA-YCV(J)*ZCV(K) 

FLOW-AREA*U(2JJC)*RHO(UJC) 

IHFF-AREA*GAM(lJJCyXDIFa) 

CALLPROFIL 

IF(ABS(FLOW)LT.  lE-20)  FIjOW=0. 
AIM(2JJC>-AOOF+AMAX1(ZEROJ1jOW) 
FLOW=AREAnj(LUJCfRHO(LUJO 
DIFF-AREA*GAM(LUJCyXDIF(Ll) 

CALLPROFIL 

AIP(L2JJK)^AOOF+AMAXl(ZEROJLOW>fLOW 
610a»4TINUE 
C  WRnE(**)*TC0EFF3  ' 

D061 11-2X2 
D061IK-2XI2 

COEFFICIENTS  ANORIH  AND  aSOUIH 
AREA-  XCV(irZlCV(K) 

FLOW-AREA*V(L2JO*RHOa.i;C) 

DIFF-AREA«GAMaiXyYDIF(2) 
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CALLPROFIL 

AJM(UK^ACOF+AMAXl(2EROfLOW) 

FLOW»AREA*VaMlJC)*RHOaMlJC) 

DIFF-AREA*GAM(LMlJCyYDIF(Ml) 

CALLPROFIL 

AJP(LM2JC)“Aa)F+AMAXl(2EROJLOWH=ljOW 
611<X»JTINUE 
C  WRrnE(**)'TOOEFF4 ' 

006121=20^ 

DO  612^2^12 

CCOnaENTS  AOUT  AND  Aihrro 
AREA*  YCV(J)*XCV(I) 
FLOW=AREA*W(UJ)*RHO(lJ.l) 
DIFF-AREA*GAM(U,  1>ZDIF(2) 

CALLPROTL 

AKM(U^)*ACOF+AMAX1(ZEROJLOW) 

FIj0W-AREA*W(U.N1)*RH0(UJ^1) 

DIFF=AREA*GAM(U>IiyZDIF(Nl) 

CALLPROIL 

AiaUJJ«)=AO(y+AMAXl(ZEROJljOWHTjOW 
612  CONTINUE 
C  WRrrE(**)'TCCEFF5  ’ 

D03987MX1 

D03987J=lJk!l 

D03987K=1J^1 

VOL*YCV(J)*XCV(I)*ZCV(K) 

APT>«HO(UXyDT 

AP(UJC)*AP(UJC>APT 

CCWJK)<XM^JC><-APTTT)(UJC) 

AP(LJ4C>-(-AP(UX)*V0L+AIP(LJJC>+AIM(UJC) 

1  +AJP(UJC>^AJM(U JC)<-AKM(UX><-AKP(IJJC)) 

2/RELAX(NF) 

apKijJcy^ijJc) 

00N(UJCHX>N(UJC)*\AX-+REL*APttJJC)*F(lJ4W 


C 

C  SAVE  OLD  VALUES  O' T  FOR  CX»JVERGENCE  CHECKING 
C 

T0LD(UK)=T(UjC) 

3987  OWnNUE 

«***•«**«**«*«****«•*•*«***««««««**«*«*«******««****««*« 

CALL7DMA 

C 


600  CONTINUE 


C  CHECKFOROONVERCENCEINTHISTIMESTCP 

C  CXX4VERGENCE  BASED  ON  CHANGE  IN  TEMPERATURE  BETWEEN  SUCCESSIVE 
C  IIERATK^S  AND  CM  OVERALL  ENERGY  BALANCE 
C 
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ICON\M) 

IFCndUEQDTHEN 

rrERi-rreRi+1 

ELSE 

C  FlhD  MAXIMUM  VALUES  FOllUSITERATK^ 

TMX-0.0 

UMX-0.0 

VMX-0.0 

WMX-0.0 

D0689Ma-l 

D0689J-1M1 

D0689K-i;^l 

TMX»AMAXl(ABS(T(UJC)).TMX) 

UMX-AMAX1(ABS(U(UJC)),UMX) 

VMX»AMAX1(ABS(V(UJC)XVMX) 

WMX-AMAX1(ABS(W(UJC)XWMX) 

689  oormNUE 
LCLTMX-0.0 
006901=2X2 
DO690J=2>42 
DO690K-2JO 

C  CALCULATE  RELATIVE  CHAN®  IN  TEMP  FROM  LAST  ITERATICW 
DELT-ABS((T(UX>-TXXD(UJC)yiMX) 
DELU=ABS((U(UJKHXXD(UX)yUMX) 
c  DELV-ABS((V(UJO-VOlI)(UJC)yVMX) 
c  DELW=ABS((W(UJC)-WOtD(UJC)yWMX) 
DELTMX»AMAX1(DELTMXJ5ELT) 
DELUMX«AMAX1(PELUMXX>ELU) 
c  IXLVMX=AMAX1(DELVMXJ[M.V) 
c  1XLWMX=AMAX1(E«LWMXJ)ELW) 

690  CXWTTNUE 

ifl;tiine.gL0.4.and.tol.gt-100)golD  1023 
IF(DELTMXGT£PST)  GOTO  691 
c  ifl;ilerl.k.20)gDto691 

C  CHECK  ENERGY  BALANCE  AFTER  ESLTMXCRTIERIA  MET 
CALLNRGBAL 
IF(ABS(EhBAL)iJE.  1 .5)  THEN 
1023  ICCWV=1 
nERL=nERl 
ITERl"^ 

ENIXF 

691  nERl=nERl+l 
IF(ABS(ENBAL).GE.200)  THEN 

C  nERATlCX^  AREIXVERCHNG:  TERMINATE  RUN 

WRTIEC*,*) '  NVERjONG  riERATTWS:  RUN  TERMWATEO 
WRnE(*,*) '  TRY  SMALLER  RELAXATTON  FACrC»S' 

STOP 

ENDIF 

IFCTIERIGT.  10000)  THEN 
CALLNRGBAL 

c  WRnE(4*)XYZXU.YVZW,UO.VO.WOJ».TO 
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WRrre(«  «) '  EXCEBXD  MAX  hO«4BEROF  FreRATK^IS  PERTIME  STEP 

WRnE(*  .•)  ■  PROCalAM^HU^IINA'^NG 

WR^E(*•)•TIME-'.TIME 

WRnE(*  •)  •  DELTMX-^  JJELTMX '  ENBAL  =‘£NBAL 
STOP 
ENDIF 
ENDIF 

CALLNRGBAL 

c  WRriE(l3*)  FIERI  J)EL^MXJ>ELlJMXXa-VMXJ>ELWMXJE^fflA^ 
IFaCXJNVEQOTOEN 

C  TIME  STEP  IS  NOT  CONVERGED;  UPDATE  BOUNDARY  CONDmC»JS  AND  ITERATE 
I)riiil*.'nME-  ’.TIME.’  ECLTMX  ‘JKLTMX 
pnit*.ENBAL;  'ENBAL'  ITERl.  'TIERI 
pfiit*,T(4^.6).T(12,7,9).T(7.13.15) 

CALLBNmY 
GOTO  99 
ENDIF 

C  TIME  STEP  IS  CC»/VERCED;INCREMEOT  TIME  AND  ITERATION  COUNTERS 
100  TTME-TIME-HDT 
riER-lTER+1 
IF1(TIMECE.TLASD 
I  LSTa>=.TRUE. 

RETURN 

end 

C 

£  »****«*«*«««*«««•»«••••*«•*«»«•**•»•««*«*««***«****•«****«***«•*»«• 

subiKutinelNCRV 

C 

LOaCAL  LSOLVEaJWOTOBLKJETOP 

Ca^1MC»I  F(35^,35,5XP(35^,35XRH0(33^,35XGAM(35^35X 

1  CONa5a5.35).AKP(35a5,35XAKM(35^35XAP(35^^5X 

2  AIP(35^,35),AIM(35^,35),AJP(33J2545),AJM(35^35) 

COMMON  deai(35^,35),delh0(35^35).epsi(35^35), 

3  X(35).XU(33)^IF(35)^CV(35)^CV^5X«LcRalaia<. 
4Y05),YV(35XYDIF(35).YCV(35).YCVS(35).tineIl,tprev(35^,35X 

5  Z(35XZW(35XZDIF(35XZCVa5)ECVS(35).ap0(35^35X 

6  YCVR(35XYCVRS(35XARX(35).ARXK35).ARXJP(35Xapip5^.35). 

7  R(35XRMN(35XSX(35XSXMN(35XXCVI(35XXCVIP(35). 

8  YCV](35XYCVJPp5XZCVK(35XZCVKP(33Xsl(35^35) 

OOMM(X4  DU(35^^5XDV(35^^5XDW(35^J5XFV(35XFVP(35X 

1  FX(35XFXM(35XFY(35XFYM(35XPT(35XQr(35XTOLD(35^,35X 

2  F2(35XFZM(35X  VHAT(35^.35XWHAT(35^,35XUOLD(35^,35) 
COMMC»^/D^X/RELAX(13XLiTUOT(13XLKJC(l  IXNTIMESCIOX 
11EOLX^10XTIMEX)T^YLZLSRHOCON7ERO,TLAST. 
2NF>BWAXJ^J4RHOJ«AMLliZL3>Il>CAOa^i;CJ13. 

3ISTJSTJOTJIERMSr, 

41FREFJFREFJCFREFMODE 

OOMMC»4/HEADINmTLE 

CHARACTER*10  TnLE(l3) 

COMMONCNIIASTOPJCALLISTOP 

COMMONAXJNViyEPSUEPSVEPSWEPST4CONVJIERl.TO(35^,35XENBAL 
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I  UO(35^35XVD(35^.33),W0(35a5J5)jreRL 
a}MMC3N«ESIDAlMAX(13)JRESID 
oaiiiin»f)icts4dbioe(35^J5)oii>ce0^.25^5X2fi^ 

COM^Oti/SOIfiaSMA^ 

CX3MM0N/a]EF/FLCWX>IFF>CX^ 

DIMENSION  U(35^35).V(3S^.35).W(35^.35)KX35^.35) 

DIMENSION  T(35^;}5) 

BQUIVALENCE(F{l.l.l,l).U(l.l.l)).(F(l.l.l^XV(l.l.l)). 

1  (F(l.l.UXW(l.l.l)X(F(l.l.l.4XPC(l.i.l)) 

2  .(F(l.l,1.5XT(l,l.l)) 

TIME  STEP  IS  CXJNVERCED;  INCTEMENT  VARIABLE  ARRAYS 
rewindCT) 

DO700I=lJ-l 
DOTOOJ-l^fl 
DOTOOK-LNl 
TO(UJC)-T(UJC) 

UO(UJC)-U(UJC) 

VD(UJO-V(UJC) 

W0(U4C)-W(UJC) 
ddh0(ijjc)-<idh(ijjc) 
l|jrev(iJjc>^ij4c) 

700  CXKnNUE 

WRTIEC?  ,*)  XYZXU.YV7W,U,V,WJ».T 
RETURN 
END 

(XCXXXXXXXXXXXXXXXXXXOCCCXXXXXXXXIXCCCCCXXXXXXXXXCCCCXXXXXXXXXXXXX: 
SUBROUTINE  UCX)F 

CCXXXXXXmXCC(XCCXXXCXXXXXXXXXXXXXXXXXXXXXX:(X(XCCCCXXXXXCCCXXCCXX: 
LOGICAL  LSOL VEJJTUNTXBLKJLSTOP 
COMMON  F(35a5.35^35^  J5XRHO(35^,35XGAM(35^.35X 

1  C(»l(3Sa5,35),AKP(35a3.35XAKM(3545^5XAP(33^,35X 

2  AIP(35a5.35XAIM(35a535XAJP(35a535XAJMC35^,35) 

COMMON  dell<33^;j5X<WW03^45Xq»i(33^;j5X 

3  X03XXU03XXDIF03XXCV03XXCVSQ3Xlkfl>3lafcnL 

4  Y(33XYV(33XYDIF(35XYCV(33XYCVS(35Xtineft,tlJrcv(33^,33X 
3  Z(33XZW(33XZDIFq3XZCV(33XZCVSO3Xap0(33^;J3X 

6  YCVR(33XYCVRS(33XARX(33XARXK33XARXJP(33)^1P3^,35X 

7  R(33XRMN(33XSX(33XSXMN(33XXCVK33XXCVIP(33X 

8  YCVK33XYCV]P(33XZCVK(33XZCVKP(33XsK33^.33) 

COMMON  DU(33^.33XDV(33^.33XDW(33^^3XFV(33)  JVP(33X 

1  R<(33XFXM03XFYO3XFYM(33XPT(33XQrO3XTC)LD(33^,33X 

2  FZ(33XFZM(35XVHAT(33a5,33XWHAT(334535XUOLD(33^,33) 
OOMMON/INDX/RELAX(13XLFNNF(13XLBLK(UXNnMES(10X 
1LSOLVE(10XTIMEJ)T^YL,ZLARHOOON7ERO.TLAST. 
2NF;«TVlAX>n»J«HOJ4GAMJLlX2X3AIl>!2>OJ41WO, 

3ISTJSTjCSrjIERJLASr. 

4IFREFJPREFJCFREFJ40DE 

0a4MCX«HEADlNmiI£ 

CHARACIER*10  TTILECD) 

COMMCX^Cm/LSTOPJCALUSTOP 
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U  O  o  U 


(XMM0N«ESa>RMAX(13)JRESD 

CX)MMO^KX)NVl/EPSUfPSV^W£PST4CONVJTERl.TO(35^5,35)£NBA^ 
1  UO(35^;}5).VO(35^.35).WO(35J5.35)4TE3U. 

(XXvlM(»I/SORCySM^ 

CX3MMON«X)EF/F1jOWJ)IFF>VCX3F 

DIMENSION  IK35^,35).V(35^.35),W(35^.35)JQ35^5.35) 

ooiiun)a^G3ioB^xfi]ne(3S^5.3S),yfoiQe(354545);z&ice(35^^ 

DIMENSKW  7(35^3,35) 

EQUIVAl£NCE(F(l,l.l,l).U(l.l.l)).(F(l.l.UXV(l.l.l)). 

1  (F(l,1.1.3XW(l,l.l)X(F(l.l,l.4)J>C(l.l.l)) 

2  .(F(1.U,3XT(1,1.1)) 


ENTRY  UCOF 

cxonciEhrrs  for  the  u  equation 

c 

CALL  reset 
NF*1 

IF{.N0T.LS0LVE(NF))  GOTO  100 

IST-3 

JST-2 

KST-2 

CALLDIFFUS 


REL=1.'RELAX(NF) 

C  WRnE(21*)'U:COEFFr 
DOl03I»34J 
DO  103  >=2>!2 
DOI03K=2>Q 

CCCFFlCIENrS  AEAST  AND  A  WEST 

FLPHja+ux)^Fxa+iri»«(i+u4C)+fXMa+irRHO(U4C)) 

FLOW-YCV(J)*ZCV(K)*0.5*(FL+fLP) 

DIFF=YCV(J)*ZCV(K)*GAM(UJCyXCV(D 

CALLPROFIL 

AIMa+UJC)=ACC»=^+AMAXl(ZEROJlJOW) 
AIP(U4C)=AIMa+UJCHTjOW 
CCEFFICIENrS  ANORTH  AND  ASOUTH 

FL=XC\^*V(I4+IJC)*(FY(3fl)*RHO(LH-lJCHfYM(J+l)*RHO(UJC)) 

FL^^XC^^-I)Wa•I4+I40*(FY(^fl)•RHOaW4+^IJC>+FifM(3^-l)• 

IRHOa-UJO) 

GM-GAM(UJC)*GAM(U+14C) 

1  /(YCV(J)*GAM(I>IJC>fYCV(J+l)*GAMaJ4C)+ 

2 1.(»B-20)*XC:VI(D 
GMM=CAMa*U4C)*GAM(I-U+lJC) 

1  /(YCV(J)*GAMa-I>IJC>^-YCV(J+ir 

2  GAMa-1  J4C>f  l£-20)*XCVIPa-l) 

MFF*ZCV(K)*2*(GM+GMM) 

FLOW-ZCV(K)*(FL+FLM) 

CALLPROTL 
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AMU+1X>-ACOF+AMAX1(ZEROJLOW) 
AJP(UX)-AJM(I^1JC>FL0W 
COSnOENrSAIN  ANDAOUr 

FlJ-XCVI(I)*W(UX+l)*(FZ(K+I)*RH0(UJC+l>4fZM(K+iri»«(UJC» 

Fii^xcMPa-i)*waMjx+i)^FZ(K+i)*iwoa-ux+i>+FZM(^ 

IRHOa-UJC)) 

CiM=<3AM(UJC)*GAM(UJC+l) 

\  /(ZCV(K)KiAM(UJC+l)+^ZCV(K+l)*ClAM(UX>^ 

2  l.(«-20)*XCVI(D 
C»4M-<jAMa-UJC)*GAMa-UJC+l) 

1  /(ZCV(K)*GAMa-UJC+l)<-ZCV(K+l)* 

2  GAMO-l  JJCK  l£-20)*XCVIPa-l) 

DffF-YCV(J)*2.*(C3M+<aiM) 

FU)W»YCV(Jr(FL+FLM) 

CALLPROFIL 

AKM(U4C+1)«ACX)F+AMAXI(ZER0JL0W) 

AKP(UJC)-AKM(U4C+l>fLOW 

103  (DONTINUE 

C  WRrrE(**)'U:CX)EFF2' 

D0104>°2M2 

DO104K-2^2 

(XenOENTS  AEAST  AND  AWEST 
AREA*YCV(J)*ZCV(K) 

FlJOW-AREA*RHO(UJCrU(2JJK) 

DIFF-AREA*GAM(UJC>XCV(2) 

CALLPROFIL 

AIM(3JJC)-A(X)F+AMAX1(ZER0J=IX)W) 

FlJOW-AREA*RHO(LUJC)n](LUX) 

DlFF«AREA*GAM(LUJCyXCV(L2) 

CALLPROFIL 

AIP(L2JJC)=ACa^AMAXl(ZEROJLOW>FLOW 

104  CX^mNUE 

C  WRnE(V)'U;CC*FF3' 

DO105I-3L2 

DO105K>2J42 

CXXFFIOEKrS  ANCXOH  AND  ASOUTH 
FL=XCVKD*V(UJCrRHO(llJX) 
FLM-XCVIPa-l)*va-lLK)*RH0a-1.140 
FLOW-ZCVaC)*(FL4fLM) 
<»4-XCVT(I)*GAMaiJC)+XCVIPa-l)*GAMa'IJ 
DIFF“ZCV(Kr<jMnrD!Fa) 

CALLFRCX^ 

AJM(ILK)=ACOF+AMAX1(ZEROJLOW) 

FL“XCVKI)*V(LM1JQ*RH0(LM1JQ 

FLM-XCVIPa-l)n^a-lA«JC)*RHOa-l>«JC) 

FLOW=ZCV(K)^-H!LM) 

GM-XCVI(I)*GAM(I>llJCHXCVIPa-l)*GAMa*^ 

EHFF*ZCV(K)*GM/YDIF(M1) 

CALLPROFIL 

AJP(I>i2JC)=ACXy+AMAXl(ZEROJLOWHLO 

105  CONTINUE 
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p 


C  WRrrE(**)’U  C0EEF4' 

DO  106 1-3X2 
DO1063-2X12 

CXJEFFICIENTS  AIN  AND  AOirr 
FL-XCVKD*W(U^)*RHO(U.l) 
ixM-xcvipa-i)*wa-u^)*RHoa-u,  1) 

FLOW=YCV(J)*(FL+FLM) 

<»^XCM(I)^jAM(U.  l>^XCVnXI-l)»GA^  1) 
DIFF=YCV(J)*C»ia)lF(2) 

CALLPROFIL 

AKM(U^)-ACOF+AMAX1(ZEROJLOW) 

FL-XCVKD*W(UJ^1)*RH0(UJ<1) 

FLM-xcvipa-irwa-u>ii)*RHoa-ijj4i) 

FLOW-YCV(J)*(FL+FLM) 

GM=XCM(I)HJAM(U^mXCMPa-l)'HjAMa-UJ^l^ 

DlFF=YCV(J)Ha4/ZDlF(Nl) 

CALLPROFIL 

AKP(U>12)=ACX^+AMAXl(ZEROiaX)W>FLOW 
106  CC»mNUE 
C  WRrrE(V)'U:CX*FF5’ 

DO107I-3JL2 

DO107J-2JM2 

DOI07K-2X12 

VOL-YCV(J)*XC:VS(I)*ZCV(K) 

AIT=(RHO(UW*XCVI(IHlHOa-U3C)*XC^ 

1/(XCVS(D*DT) 

AP(UJC)-AP(UX>APT 

CX)N(UJCHX)N(UX>^■AP^•lK)(UJC) 

AP(UJC)- 

1  (-AP(UJCrVOL+AIP(UJQ<'AIM(UJC)<-AJP(UJO+AJM(UX) 

2  +AKP(UJCKAK\VJJC)) 

3/RELAX(NF) 

mm*****t‘**»***»********»*******************************souifX\am 

0m*m*mm*******»******m**********m»*************»***************» 

CX>N(U^)-CX)N(UJCrVOL+REL*AP(UJQ*lWJJC) 

DU(UJC)=VC*-«DIF(D 

DU(UJCH>U(UJCyAP(UJC) 

i07a»mNUE 

IFaCALLEQ.l)THEN 

C  TEMPCaLWY  USE  (yPCOJ)  TO  STORE  UHAT _ 

D0151K=2J« 

DO  151  >2JV12 
DO  151 1=3X2 

PC(UXHAIP(UJC)*Ua+UJC)+AIM(UJC)*Ua-UJC) 

1  +AJP(UJC)*U(U+lJC)+-AJM(UJCrU(U-lJC) 

2  +AKP(UJC)*U(UJC+1HAKM(UJC)*U(UX-1) 

3  4CON(UJC)yAP(U4C) 

151  ooainue 

endif 

100  C(»mNUE 
RETURN 
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end 


C 

Q  •«««**••«*«**««««««««*««•«*•««••«•«•«•«•«*«•••«••*«*••«««•««•**«•»*« 

aimutineVGCX^ 

LXX3C AL  LSOL VE^ilUKrjLBUaSTOP 

COMMON  F03J3;35.5)«35^.35)JlHO(35^35XGAM(35a5,35X 

1  CON(35a5^5).AKP(35^,35).AKM(35^,35)^5^.35X 

2  AIP(35^^5),AIM0545,35).AJP(35^35).AJM(35^.35) 

COMMON  (Uh(35^.35).ddhO(35^,35),«|JSi(35^,35), 

3  X05)WJ3X®IF(35XXCV(35)^CVS(35X&.cp.alaienl. 
4Y(35).YV(35),YDIF(35XYCV(35).YCVS(35Xtinelt,lpiev(35a5.35). 

5  Z(35)7W(35)7DIFO5)^(35)7CVS(35)4^j0(35J5.35X 

6  YCVR(35XYCVRS(35XARX(35),ARXK35).ARXIP(35Xapl(35^.35X 

7  R(35XRMN(35XSX(35XSXMN(35XXCVI(35XXCV1P(35X 

8  YC\^35XYCV]P(35),ZCVK(35)7CVKP(35)^35^.35) 

COMMON  DU(35^,35XDV(35^^5XDW(35^35XFV(35XFVP(35X 

1  FX(35)J=XM(35)JfY(33)JYM(35)JT(33XQTa5XTOLD(35^,35X 

2  FZ(35XFZM(35XVHAT(3543,35XWHAT(35^35XUOLD(35^^5) 
COMMON/INDXdtELAX(13Ui'IUNT(13),LBLX(llXNnMES(10X 
lLSOLVE(10XTIME;DT;a-.YUL.SJtt«CON7ERD.TLAST. 
2NFJ^FMAX>IP.NRHOJ^GAMJ-li2JJJM!A!2AC;^lW 
3ISTJSTJCSTJTCRMST. 

4IPREFJPREFJCPREFMODE 

COVlMC^j/HEADINrrniE 

CHARACreR*10  TnLE(13) 

COMMONChniASTOPJCALUSTOP 

CC»4MC»<m4VL®»SU,EPSVJEPSW£PST.ICC»4VJIERl.T0(35^535)£^^ 
1  00(3505.35), VO(35^,35XWO(35^,35XnERL 

C0MM0NA(ESID/RMAX(13)JRESID 
CX)MMON/SOR(:/SMAX,SSUM 

ooinnK)iidbf0^x&)foe(350505)^lbfee(3S0505),2£xoe(3S05,3S) 

COMMONADOEF/FLOWJDIFF^^COF 

DIMENSION  U(350305X  V(3505,33XW(3503.35XPC(3503.35) 

DIMENSION  7(3505,35) 

EQUIVALENCE(F(l,i,l,lXU(l.l.l)),(F(I,I,10XV(l,l,l)X 

1  (F(l,l,1.3XW(l,l.l)X(F(l,l,1.4)J>qi,l,l)) 

2  .(F(1.1,1,5XT(1.1,1)) 

C 

CXJEhFKJltNTS  PC»THE  V  EQUATION - 

C 

CALL  RESET 
NF=2 

IF(.NOTLSOLVE(NF))  00  TO  200 

IST»2 

JST*3 

KST*2 

CALLDIFFUS 

REL«1.-RELAX(NF) 

C  WRnE(V)'VOOEFFl' 

DO203>3>C 

DO203I*2iO 
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DO203K-2JC 

CX)EFFraENrS  ANORTH  AND  ASOinH 

FL-V(UJK)^(JrRHO(UJC>tfYM(J)*RHO(^ 

FU^V(U+lX)^(J+l)*«HCXI>lXHTM(J+l)*RHOaJX^^ 

FU)W-XCV(I)*ZCV(Kr0.3*(FL+FLP) 

DIFF^XCV(irZCV(K)*GAM(UJCyYCV(J) 

CALLPROFIL 

AJM(I^1JC>-ACOF+AMAX1(ZEROJ1jO>^') 
AJP(UJC)-AJM(U+lJChFLOW 
(XEFFKaENTS  AEAST  AND  AWEST 

FL-Y<:vKJrua+ujonFxa+i)*RHoa+ijx)+FXMa+irRH^ 

FlJ^YC\aP(J-lW+l>lJK)*(FXa+l)*RHOa+14-JX)+fXMa+l)* 

1  RH0(U-1  JC)) 

avWiAM(UJC)*GAMa+UJC) 

1  /(XCV(I)*GAMa+iJJ^>^X(:^'a+l)*GAM(UJC>f 

2  l.0E-20rYCVJ(J) 

1  /(XCV(D*GAM(I+U-lJC>fXCVa+l)* 

2  GAM(U-lX)<-l£-20rYCVIP(J-l) 

DIFF-ZCV0C)*2.*(GM+<»1M) 

FL0W-ZCV(K)*(FL+FLM) 

CALLPROFIL 

AlMa+lJ4C>‘ACC»f+AMAXl(ZEROJLOW) 
AIP(UJC>«AIMa+lJ4C>FLOW 
OMFFICIENTSAIN  AND  AOUT 

Il-YCVXJ)*W(U4i+l)*(FZ(K+l)*RHCKUJC+l>+FZM(K+irRHOaJ4C)) 

FLM-YCV]P(J4)*W(U-lj;+ir(FZ(K+l)*l«CKl>lJC+lH=^^ 

1RH0(U-IJC)) 

GMKjAM(UJC)*GAM(UJC+1) 

1  y<ZCV(K)*GAM(UJC+l>fZCV(K+irGAM(UJCK 

2  l.(E-20)*YCV](J) 

GMM><3AM(U*lJCrGAM(U-lJC+l) 

1  /(ZCV(KrGAM(U-l4C+l)+^ZCV(K+l)* 

2  GAM(IJ-14C>f  1£-20)*YCVJP(;-1) 

DIFF=XCV(I)*2.*(C^!+<jMM) 

FLOW=XCV(I)*(FL+FLM) 

CALLPROFIL 

AKM(IJJC+1^AOW^+AMAX1(ZEROJLOW) 
AKP(U4C)“AKM(UJC+l>FLOW 
203  CONITNUE 
C  WRrre(**)’V;COEFF2' 

D0204I»2X2 

DO204K»2^ 

0(£FnCIENrS  ANORIH  AND  ASOUTH 
AREA«XCV(D*ZCV(IO 
FLOW^AREA*RHOaiJCrV(L2JC) 

DIFF*AREA*GAM(LlJCyYCV(2) 

CALLPROFIL 

AJM(L34C)»ACOF4-AMAX1(ZEROJLOW) 

FljOW»AREA*RHOaMljC)*V(LMlJC) 

DIFF=AREA*GAM(LMlJCyYCV(M2) 
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CALLPROFIL 

AJP(IJ»42JQ-ACC)F+AMAXUZEB0J1DW)^ 

204  OONTINUE 

C  WRrrE(**)’V  COEFF3' 

DO20S^3>Q 

D020SK-2JC 

CCCFFIOB^  AEAST  AND  AWEST 
FL-YCYJ(J)*U(2JJC)*RHO(UJC) 
FU4-Ycvjp(J-iruaJ-ixrRHO(i>ijc) 
FL0W-ZCV(K)*(FL+FLM) 
GM-YC\^J)KiAM(lj;C)fY(:V]I^>l)*GAK^ 
DIFF-ZCV(K)*GN«XDIF(2) 

CALLPROFIL 

A1M(24JC)-ACOF+AMAX1(ZEROJIjOW) 

FL-YCVJ(J)*U(L1JJC)*RH0(LUJ0 

FLM-YCVJP(J-l)nJ(Ll>lJC)*RHC)(LU-!JC) 

FLOW-Zt:V(K)^+FLM) 

GM-YCV3(J)HiAM(LlJX>^-YCVIP(J-l)*GAM(LI>lX) 

DIFF-ZCV(K)*®4/XDIF(L1) 

CALLPROFIL 

AIP(L2JJC>-ACXy+AMAXl(ffiROJlX)W>FLOW 

205  CONTINUE 

C  WRnE(*.*)'V:C0EFF4' 

DO  206 1-2X2 
DO  206  >3X0 

COEFFICIENTS  AIN  AND  AOUT 
FL-YCVJ(J)*W(UXrRHO(U,l) 
FLM-YCVJP(3-lfW(U-U)n<H0(U-Ll) 
FLOW-XCV(D^4FLM) 
Ca>YCNa(J)KiAMajJ>+-YCNai\J-l)*GAM(lJ-l,l) 
DIFF-XCV(D*GN«n)IF(2) 

CALLPROFIL 

AKM(UX>-ACOF+AMAX1(ZEROJ=LOW) 

FL-YCVJ(J)*W(UJJ1)*RH0(UJ41) 

FLM-YCVJP(J-irW(lJ-l^l)*RHO(U-l^l) 

FLOW-XCV(D*(pL-HXM) 

(aw>YCVKJ)*GAl^J4l)fYCVJP(J4)*GAM(LXl^ 

I»FF-XCV(irGM2DIF(Nl) 

CALLPROFIL 

AiaWJC)-ACOF+AMAXl(ZEROFLOW>«jOW 

206  CONTINUE 

C  WRnE(**)'V:O0EFF5' 

DO  207  >2X2 
DO  207  >3X0 
DO207K-2XI2 

VOX-XCV(D^CVS(J)*ZCV(K) 

APTMRH0(UJC)*YCV3(J>+RH0(LJ-lJCrYCVIP(Xl)) 

y(YCVS(J)TO 

AP(U40AP(UJK>APT 

CON(UJC)-CON(UX>AFT*VO(UX) 

AP(UJ0- 
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n  o 


r 


1  (-AIWJC)*V0L+AIP(UJC>fAIM(U4C>+AJIVJJC)fAJM(UJC) 

2  +AKP(UJC)f AKM(U40) 

3flffiLAX(NF) 


»♦»♦»»»»»»*  W»««««**M»**«*»*«**«****«***«W**»**»*»*ap.nrr^inT?« 


CX)WX)<X>N(UJC)*V(DL+REL*AP(UJC)*V(UX) 

DV(UaC)F=VOI/miF(J) 

DV(UJC)>DV(UJCyAP(UJC) 

207CXKnNUE 

IFaCALLBQDTHEN 

0080991-20^ 

D08099K-20J2 

DO8099J-3A12 

8099  VHAT(UJCMAIP(UJC)*Va+UJC)fAIM(UjK)*Va-104C) 

1  +AJP(UJQ*V(U+IX)+-AJM(UJC)*V(U-1X) 

2  +AKP(U0C)*V(UJC+1>4-AKK^UJK)*V(U0:*I) 

3  ■KX)N(UX)yAP(UJC) 

ENIMF 

200(X»<mNUE 
RETURN 
end 

M****************************************************************** 

sutmudneWCXIF 

UX3CAL 1^ VEOJ>RlNrOBLKO.STOP 
COMMON  Fa3,25,35,5XP(3305.35)J«!0(334535).GAM(35^.35X 

1  CX)NO3^33).AKP(33O3O3)0Ua^3O343)UHPO3O^^^ 

2  AIPO3a333)>UM05O3.33)OVJPO3O333)>UM(33O5.33^ 

COMMON  delh(33^.33Xdeili0(33^43),epsi(33^33). 

3  X(33)XUO3)0a3IFq3pCCVa3)0(C\^3Xtk,cp,alato^ 

4  Y(33XYV(33XYDlF(33XYCV(33XYCVS(33).tindt.lprcv(3303,33X 
3  Z(33);W(33)0I3IFO3)^(33).ZCVS(33Xa|j0O3^;J3). 

6  YCVR(33),YCVRS(33XARX(33),ARXJ(33XARXJP(35Xapl(3503,35X 

7  RO3)«y4NO3)^35XSXMbK33)0CCVI(33XXCVIP(33X 

8  YCV103XYCVlPa3),ZCV1C<33).2CVKIV3)^^3 
COMMON  DU(3303.33)JDVa303.33)J>W(3303^3XFV(35XFVP(33), 

1  FXO3)J»XM03)JTO3)EYMp3)JT(33XQrO3XTOLDO3O333^ 

2  FZ(33XEa4(33),VHAT(33^^3).WHAT(3303 J3XU0lD(33O3^3) 
OOMMON/I^OX/RE1^13)OHU^^'(13)0£LK(1  lXNnMES(10X 
11^VE(10XT1MEJDTOCUYL2US1«OCON7ERO.’ILAST, 
2NF04FMAX0ffJflW0iGAMOlX20J041JV!2J^OI10i^ 
3ISTOSTJCSrjreROAST, 

4IPREFJPREFJCPREPM0DE 
OOMMON«EADO^mLE 
CHARACIER*10  TTn£(13) 

OOMMONCNrULSTOP4CALL4STOP 

C0MM0Nm4Vl/EPSUJEPSV^W£PST4C0NVJIERl,TD(3303,33XENBAL. 
I  UO(33^^3XVO(330343XWO(330333)JIERL 
OOMMON/RESID«MAX(13XIR£SID 
OOMMON/SORC/SMAX,SSUM 
COMMON/COEF/F1jOW43IFFOVCOF 
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DIMENSION  U(35  J5.35XV(35a535XW(35a535)JPC(35^;35) 
ooiiimaa^bRX^faa03.2S^S),yface(35^^SX2fcra(3S^.35) 
DIMENSION  T(33^,33)  * 
BQUIVALENCE(F(1,1.1,1XU(I.1.1)X(F(1.1.UXV(1.1.1)X 

1  (F(l.l.l.3XW(I.l.l)X(F(l.l.l,4XPC(l.ia)) 

2  ,(Fa.l.l^T(l.l.l)) 

C 

CX3ErFICIENrSPOR™EWEQUA7X)N - 

C 

CALL  RESET 
NF-3 

IF(.NOTLSOLVE(NF))  GO  TO  300 

isr-2 

JST-2 

KST«3 

CALLDIFFUS 

REL-1.-RELAX(NF) 

C  WRIIE(**)'W;COEFFr 
DO303K=3>Q 
DO3031-2M2 
DO  303 1-2X2 

COEFFICIENTS  AIN  ANDAOUT 

FL-W(UJC)*(FZ(K)*IM0(UJC>H7M(K)*RH0(UJC-^ 

FLP-W(UJC+1)*(FZ(K+I)*RH0(UJC+1>+FZM(K+1)*RH0(UJC)) 

FLOW-YCV(J)*XCV(I)*0.5*(FL+FLP) 

DIFF»YCV(J)*XCV(I)KiAM(UJCyZCVaC) 

CALLFROFIL 

Aia4(UX+l^AOOF+AMAXl(ZEROJLOW) 
AKP(UX)-AKM(UX-*'1>FL0W 
COEFFICIENTS  ANORIH  AND  ASOUm 

FL-ZCVK(KrV(LJ4-lX)^(J^-I)*RHO(LF^lX)HTrM(3+irRHO(LJX)) 

FLM-ZC^XP(K-lrV(L^^IX-l)^FY(^^■l)*RHO(L^lX-i)+FYM(ffI)• 

1RH0(UX-1)) 

CiM-GAM(UX)*GAM(U+IX) 

1  /(YCV(J)*GAM(IX^1X>+-YCV(J+I)*GAM(LJX)<- 
2 1.(E-20)*ZCVKaC) 

GMNMjAM(UX’l)*CiAM(l^lX-l) 

1  /(YCV(J)*GAM(LH-IX-l)+^YCV(J+ir 

2  GAM(UX-1>^1£-20)*ZCVKF(K-1) 

DIFF-XCV(I)*2.*(GM+GMM) 

FL0W-XCV(D*(FL4fLM) 

CALLFROFIL 

AJMa>lX)-AOOF+AMAXl(ZEROJLOW) 
AJF(LJX>-AJM(U+IX>FL0W 
COEFFIOENrSAEAST  AND  AWEST 

FL-2C\lCaC)nJ(I+UX)*(FXa+I)*RH0(I+UX>+FXM(l+irRH0^ 

FLM-zc\«F(K-i)*ua+ux-i)*(FX(i+i)*RHoa+ux*i>+fXMa+ir 

1RH0(UX-1)) 

gmhjam(ux)*gam(i+ux) 

1  /(xcv(D*GAMa+ux>+^xcva+irGAM(ux>^ 

2  1.0E-20)*ZCVK(K) 
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GMMH3AM(UJC-irGAMa-t-lJX-l) 

1  /(XCV(D*GAM(KiJJC-l>fXCVfl+l)* 

2  GAM(UJC-l)^l£-20y*ZCVKP(K-l) 
DIFF-YCV(Jr2.*(GM+CMM) 
FLOW^YCV(J)^+FLM) 

CALLFROFIL 

AIMO+UJO-ACXDF+AMAXICZEROJLOW) 

AIP(U4C)-AIMa+UJC>FLOW 

303  CONTINUE 

C  WRnE(**)W:C0EFF2' 

DO304J-2>42 

003041-2X2 

COEFnOENrSAIN  ANDAOUT 
AREA-YCV(J)*XCV(I) 
FLOW-AREA*RHO(U,l)*W(U^) 
D1FF-AREA*GAM(U.  l>ZCVa) 

CALLPROIL 

AKM(U.3)-ACOF+AMAX1(ZEROXLOW) 

FLOW-AREA*RHO(UXJ1)*W(UJ41) 

I»FF-AREA*GAM(UX1VZCV(N2) 

CALLFR0F1L 

AiaUJJC)-ACW+AMAXl(ZEIWD4TjOWHljOW 

304  CONTINUE 

C  WRnE(*.*)'W:COEFF3' 

003051-^ 

DO305K-3J« 

CCeFFIOENTS  ANORIH  AND  ASOUTH 
FL-ZCVK(ICrV(IXK>»1W0aiX) 
FlAf-2CN«P(K-irV(UK-irRH0aiJC-l) 
FLOW-XCV(D*(FL+FLM) 

C»4-ZC\«(K)*GAMaiJC)fZCVia\K-l)^JAMai4C^^ 

DIFF-XCV(D*GNVYDIFa) 

CALLPROFIL 

AJM(i;UO-AOOB-AMAXl(ZERCUXOW) 

FL-ZCVK(lCrV(l>llJO*RHO(LMlJC) 

FTA«'<ZCVKI\K-l)*V(I>flJC-l)*RHO(^ 

FLOW-XCVOynFL+fLM) 

GM-ZCVK<K)*GAMaX414CKZCVKP(K-l)HiAM(^ 

DnFF>-XCVa)*C»(VYIMF(Ml) 

CALLPROFIL 

AJP(LM2JO-ACOFfAMAXl(ZEROXLOW>FLOW 

305  CONTINUE 

C  WRnE(**)'W:C0EFF4' 

DO306K-3JO 

DO306J-2XQ 

OOEFFiaENrS  AEAST  AND  A  WEST 
aFZCVK(K)*UaJJC)*RHO(UJC) 
FLM-ZCVKP(K-irU(2JJC*l)*RH0(UJC-I) 
FLOW»YCV(J)*(FL-H!LM) 

GM-ZCVK(K)*GAM(UJC)fZCVKP(K-l)*GAM(U4C-l) 

DlFI®-YCV(JfGMXIMF(2) 
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CALLPROFIL 

A1!^JC)-ACOF+AMAX1(ZEROJLOW) 

I^RZC\«(K)nKLlJX)*RHO(LlJ4C) 

aM-ZC\«P(K.l)nj(LlJJC-l)*RHO(LUjC*l) 

FLCW-YCV(J)^+fLM) 

GA«CVK(K)*GAM(LUX)^2C>«W-l)KiAK^ 

DIFP-YCV(J)*C»4«DIF(L1) 

CALLPROFIL 

AIP(L2JX)-A<X3F+AMAX1(2EROJ1jOW>FLOW 
306  OOFTITNUE 
C  WRrrE<**)’W:00EFF3’ 

DO3071-2JL2 

DO307J-ZM2 

DO307K-3JSQ 

VOU«YCV(J)*ZCVS(K)*XCV(D 

APT-(RHO(UJC)*ZCVK(K)HWO(UX-l)*ZCVia»(K^ 

lAZCVS(K>n)D 

AP(UJC)-AP(UJC)-APT 

(MN(I44C)<X)N(UJC>fAPT*W0(U4C) 

AP(UJO“ 

1  (.AP(UJ0*VOL+AIP(UJC>+AIMaJJO^AJP(U4C>fAJ^ 

2  +AKPI(UX>^AKKKUX)) 

3/RELAX(NF) 

*«»**«**»««**»««****«**«»*«**<»*»»»**«»*«*»»*«»**»*«»*»*gouroe  tenn 
•**«*««**«•*«*««*««*««•«•••««••««««*«***•*•**•*••*••****«•**«•*« 

COI^X)HX)N(UJC)*N«-+REL*AP(UJC)*W(UW 

DW(UJC)-V0L/ZD1F(K) 

DW(UJKH>W(UJKyAP(UJC) 

307CONnNUE 

IFOCALLEQDTHEN 

DO9099t-2JJ 

DO9099J-2AC 

DO9099K-3JC 

9099  WHAT(UJCMAIP(UJ0*Wa+lJJC>vAIM(LJJ0^a-1440 

1  +AJP(UJCrW(I^lJC)4-AJM(LJX)*W(U-lJO 

2  +AKP(U4CrW(UJC+l)<-AKM(UJCrW(UJC-l) 

3  -KMNCUJOyAPdJaO 
ENDIF 

300  OOKHNUE 
REIORN 
end 
C 

stdxtntine  POOF 

LOGICAL  LSOLVE2J>Rl^^JLBLIaSTOP 

COMMON  F(354535 J)JV5^^5XRH0(35^35XGAM(35  J535). 

1  CON(35^35XAKP(35^45XAKM(35a535XAP(35^35). 

2  AIP(35a535XAIM(35^^5).AIP(33a3,35).AJM05^35) 

COMMON  deO<3S^.3S),deihO(3SaS^5Xepii(3Sa5,3SX 

3  X0S)XU05)^IFaS);CCV05)^CVSqSXlk,cp^^ 

4  Y(35).YV(35),YDIF(33).YCV(35),YCVS(35XttiidUprcv(35^^5X 
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5  ZO5XZWO5);5DFO5)ZCVO5).ZCN^5)j|j0(35^^ 

6  YCN^5XYCN^5).ARX(35).ARXJ(35),ARXJP(35)4|)1(35  J5.35). 

7  R(35);»4NP5)^5).SXMN(35);CCV1(33)^CVIP(3^^ 

8  YCVJ(35XYCVJP(35)7CVK05).ZC\W5)^35^^^ 

CXJMMON  DiJ(35a5.35XDV(35a3.35)J3W(33^^5U^'(35XFVP(35X 

1  I^5)JT(M(35)JFY05)J'YM(35)JT(35).QT05XT0IJX35^^^ 

2  FZ(35XFZM(35).VHAT(35^,33).WHAT(33^.35XUC)LD(33^5,35) 
CX3MMON/D^XA<ElJO((13)JLFRl^r^(13)JL^^ 
ILSOLVt(10XTTMEJ)T^YlASJWOOON7ERD.TLAST. 

3ISTJSTJCCTJTCRMST. 

4IFfiEFJPREFja>REFMODE 
COMMON^iEADIhmiLE 
CHARACreR*  iO  TnLE(13) 

(XJMMONOrriASTOPjCAIIJ^^ 

a)MMC»4/CC»4ViyEPfaj^PSVfPSWfIOT4CCWV.nmi.T0O5^.35)3J^ 
1  \m5;25^5),\f)OS:25^S)yM^5:25,35)JTm. 

COMMOItWES[D/mAX(l2)JB^ 

(XMMON/SORC/SMAXSSUM 

OOMMONCOEF/FLOWJHFF^OOF 

ccmmoii«bice^feiOB(3^^ 

DIMENaON  U(35^.35),Va5^.35XW(35^45XPC(35^35) 
E»MENSIC»4T(3S;ZS^S) 

EQUIVALENCE(F(1.1.1,1XU(1.1.1)X(F(1.1,UXV(1,1.I)X 

1  (F(l.l.l^XW(l,l,l)X(F(l.l,1.4XPC(l.l.l)) 

2  ,(F(1,1.1AT(1,1.1)) 

C 

<MEFFraE^^^  FOR  TOE  PRESSURE  EQUATION _ 

C 

NF-NP 

IST^2 

JST-2 

KST«2 

c  WRnE(V)'PC:COEFFl' 

DOS01J-2M2 

DO501K=2442 

AIM(2JJC)-<).0 

AIP(L2JJCH).0 

CON(244C)-RHO(lJ4QnK2J4CrYCV(J)*^ 

CON(LlJjC)iO. 

501  CONTINUE 
DO  502 

DOSQ2K-2442 

AJM(UKH).0 

AJP(I>I2JCH).0 

CON(U4C)-RHOaiJC)n'(I^)*XCV(I)*ZCV(K^ 

CON(LM1JC)-0. 

502  CONTINUE 
DO  503 1*242 
DO  503  4-24^ 

AKM(U2H0 
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AKP(UJC)-«.0 

cx)W^Hwo(uarw(u^rxcv(i)*Ycv(^ 

503  CXDKTINUE 
C 

C  WRrre(**)’PC:CX)EFFy 
DOS04K<N2 
DO504>«2J^ 

005041-202 

AREA-YCV(J)*ZCV(K) 

ARHC>-AREA*(FXa+l)*«HO(I+lJX)+FXMa+l)*RHCW 

FLDW-ARHO*PCa+lJ0C) 

IF(IEQ12)  FLOW-ARHO*U(L10JO 
IF(ABS(FL0W)1,T.  lE-20)  FLOW-0. 

AIP(UOC)^ARHCWXJa+lJOC) 

AIMa+UJC>-AIP(UJC) 

CX)N(UJC)KX>N(U0C>FLOW 

coNa+ioocHWi+uocHijow 

c 

AREA-XCV(I)*23CV(K) 

ARHC>«AREA*(Fy(J<-l)*RHO(I>lJCH?YM(W)*RH^ 

FLOW-ARHO*VHATah-lJC) 

IF(J£Q^Q)FLOW-ARIiO*VaM10C) 

IF(ABSO^W)LT.1E-20)  FLOW-0. 

AJP(UJK)-ARHOn)V(U+lJC) 

AJM(U+lJC)-AJP(UjC) 

OON(UXKX)N(UJC>«X)W 

CON(l>lJC)«OON(I>iOC>tFLOW 

AREA-XCV(I)^CV(J) 

ARHO-AREA*(F^+l)*RHO(UJC+l>+fZM0C+l)*W^ 

FLOW-ARHO*WHAT(UJC+1) 

IF(KE1Q.N2)FLOW-ARHO*W(U011) 

AKP(UOO-ARHO*DW(UJC+1) 

AKM(U0C+1)-AKP(U0C) 

IF(ABS(FLOW)XT.  1&20)  FLOW-0. 

CX3I^JQ-OON(UOC>fLOW 
CONCUOC+IHOjOW 
C  WRnE(**)'3’iOJC 

AP(UOC)-AIP(UJC>+-AIM(UJC>fAJP(UOC><-AJM(UJC) 

1  +AKP(UJC>^AKM(UOC) 

504  OCKITNUE 
RETURN 
END 

CCCCCCCCCCCCXXXXCCCCCCCXXXXCCCCXXCCCCC^ 

SUBUXJTINE  umedi 

(xxxmxcmxxxxxccxxxxxxxxxx(xccccccxxxxcxxxxx(xcmxxxcc 

LOGICAL  LSOLVEOi>IUOTO£LiaSTOP 

COMMON  F(3525254XP(352525XRHO(35a525XGAM(3525.35X 

1  C0N(3525.35),AKP052525),AKM(352535XAP(352535). 

2  AIPP52525)>UM052523)AJPO52525)>UM(3525.35) 

COMMON  delh(352525Xdelh0(352523).cpsi(3523O5). 

3  X(35)0a;O5)0a)IFO5)0(CVO5)0CCVS(350tk.cMaen^ 
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4  Y(35XW(35XYDIF(35XYCV05XYCVS05)4inelt.l|jrev0 

5  Z(35)ZW05).Z[)IFq5)^a5),ZC^rSa5);^^ 

6  YCVT<(35XYCVRS(35),ARX(35XARXJ(33).ARXIP(35Xapl(35^.35), 

7  Ra5)ja^5)^5XSXMN(35);<C\^5);CCVIIX35). 

8  YC\W5XYC^^5).ZC^/K(35),ZCViaV5)A35^^ 

COMMON  DU(35^.35XDV(3543.35)J)W(33^43)J=V(35)JVP(35X 

1  I^3XFXM(35)JYa5)JYM05)Pr(35XQrp5XTOLD(35^,35X 

2  FZ(35XFZM(33),VHAT(35^.35).WHAT(35^^5).UOLD(35a5,35) 
aDMMON/INDXmELAX(13UJWhnX13UBIJC(ll)J^nMES(10X 
ILSOLV^IOXTTMEJJT^YL^SJRHOOW.ZERO.'ILAST. 
2NFJ4r^l4X;^J«HaNGAMXli2aJJ^lM2J^J^ 

3isTjsrjcCTjiHaAsr. 

4IPREFJPRB^JCPREFMOOE 

OOMMON/HEADINnmE 

CHARACIER'IO  TniE(13) 

COMMONC^mA5IOPJCALUSTa> 

COMMON«:ONVl/EPSU,EPSVfPSW£PST.ICONVjreRl.TO(35^.35XENBAL. 

I  l»(33^,35),VO(35^,35).WO(33^J5)jreRL 

(X)MMC»4/SC»£/SMAX.SSUM 

(MMMON/COEF/FLOWJDIFFACOF 

DIMENSION  U(35a5,35),V(35a545).W(35a535XPC(35;i535) 

oonimaa«bro0SdbicB(3543.35)j«)ioe(35^;}5^ 

DIMENSia4T(35a5J5) 

EQUIVALENCE(F(I.1.1,IXU(1,1.1)X(F(1,1.UXV{1.1.1)X 

1  (F(M.1.3XW(l.l.l)X(F(l,l,l.4XPC(I,l.l)) 

2  .(F(I,I.14XT(I.M)) 

C - 

10  P0RMAT(26(lH*X3X^103X26(lJn) 

20  FC»MAT(1X4H  I  >  %.6I9) 

30F0RMAT(1X1HD 
40  FC»MAT(1X42,3X,IP.7E9.2) 

SOPORMATCIH) 

51  F0RMAT(1X,1  ='aX,7(I4.5X)) 

52  FORMATCIX.'X  ^.1P,7E9.2) 

53  F0RMAT(TH<1P,7E9J) 

54  FORMAT(lX.'J  ^JX.7(I4^) 

55  PORMAT(IX.Y  “',1P.7E9J) 

56  F(»MAT(1X.K  -’JX.7^4X)) 

57  PC»MAT(1X.Z-',1P,7E9.2) 

59  F0RMAT(1X,X  -^,2X^4) 

0*««**«*«***«***«**«***«*«****««***«*««*«*«*«*«* 

c  ENmYUMESH 

XU(2>-0. 

DX=XI;FL0AT(L1-2) 

DOlWXl 
1  XU(I)»XUa-l>+DX 
YV(2>>0. 

DY=YI;FIjOAT(M1-2) 

D02>3M1 

2YV(J>-YV(J-1>+DY 
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Twnw) 

DZ-ZW=L0AT(Nl-2) 

D03K-3JJ1 

3ZW(K)-ZW(K-1H)Z 

RETURN 

end 

(y^mm****************************************** 

subraudne  PRINT 

^«*«««««»«««**«««««««*«««««««««*«««**«*««**««**« 

LOCaCAL  LSOLVEJJWNTXBLKl^^ 

COMMON  F(35^^3^J>(33^.35XRHO(35  J5;}5).GAM(35^.35X 

1  CON(35^.35).AKP(35^J5).AKM(35a5.35XAP(35^35). 

2  An^33^.35)>UM(35a5.35XAJP(354S.35).AJM(35^.35) 

COMMON  deih(3S^,3SXdelta0(3S^.35Xq]^3S43,3SX 

3  X(35)W35)^IF(35);CCV(33)^CVS(35Xflc.C|Wdal^ 

4  Y(35XYV(35).YDlF(33XVCV(35XYCVS(35Xlinefc.lprev(35^.35X 

5  ZO5)7W(35)^IF(35)^(35XZCVS(35Xa|>0(35^35X 

6  YCVR(35).YCVRS(35XARX(35)ARXK35XARXJP(35Xapl(35^.35). 

7  R(35)JIMN(35XSX(35XSXMN(35)^CM(35);CCVIP(35). 

8  YCV3(33XYCVJI\35).ZCVTC(33),ZCVKP(35Xfl(35^.35) 

COMMON  DU(334535XDV(35^,35XDWP5^35XFV(35XFVP(35X 

1  I^5)JXM(33)JT05)JFYM05)JT(33XQm5XT0UV 

2  F2(35XF2M(35).VHAT(35^,35),WHAT(35;253SXUOLD(35^.35) 
COMMCMdNDX/RELAX(13).LPRlNT(13U^i^IlXNnMES(10). 
1LSOLVE(10XTIMEJ3T^Y1^SRHOOQN,ZERD.TLAST. 
2NFNFMAXNPJ4RHOJ4GAMJLia^J41J«42M3J4ia42J43, 
3ISTJSTJCSTJIERaAST. 

4IFREFJFREFJCFREFMQDE 

COMMON«EADINmnLE 

CHARACIER*10  TnLE(13) 

COMMON/CNllASrOPJCALLJSTOP 

COMMONm4Vl/EPSU£PSV£PSW£PST4CONV,nERl.TO(35.25.35XENBAL 

1  UO(35^.35XVD(35^35XW0(35^,35XnERL 

OOMMON^OROSMAX,SSUM 

COMMC»WX)EF/FLOWJ)IFF^COF 

DIMENSION  U(35a5.35XV(35^.35XW(35a535XPC(35^;35) 

oaininoi^bice^xfcice05J3.35X}'fiace(3SJSJSX2fiiioe(35^,3S) 

DIMENSION  T(35^45) 

EqUIVALENCE(F{I.I.I,IXU(l.M)X(F(I.l,UXV(l,l,l)X 

1  (F(1.I.UXW(I,I.I)X(F(I,1.I,4XPC(1.1,1)) 

2  .(F(l.l,14XT(l,I.l)) 


10  PORMAT(26(IH*X3XAiO,3X26(lH*)) 
20  F0RMAT(1X4H  I  >46,619) 
30FORMAT(1X1HD 
40  P0RMAT{1X424X  1P,7E9.2) 
50FORMAT(1H) 

51  FORMATdXl  ‘~^JK7(I4JX)) 

52  FORMATdXX  <1P,7E9.2) 

53  FORMATCIH  <1P.7E9.2) 

54  F0RMAT(1X*J  -■,2X7(I4.5X)) 
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55  FORMATXlX.'Y  -^.1P,7E9.2) 

56  FORMATdXK  -'^7a4,5X)) 

57  F0RMAT(1X.Z  <IP,7E9.2) 

59  F0RMAT(1X.K  =’;iXJ4) 

IF(.NOTU>RINr(3))  GO  TO  80 
gOCXWTINUE 
C 

IF(.NOTIJRINT(NP))  GO  TO  90 
C 

CONSTRUCT  BOUNDARY  PRESSURES  BY  EXTRAPOLATKW 
C 

D091K=2J0 

D091J-2AI2 

P(UJC)p<p(2JJiC)*XCVS(3>P(3JJCrXDIF(2)yXDIF{3) 

91  P(LUJCHP(l^W*XCVS(U>P(L3JJC)*XDIF(Ll)VXDIF(L2) 
D092K-2J42 

DO  92 1=2X2 

W1JC)'^(P(UK)*YCVS(3HV.3JC)*YD1F(2))W1F(3) 

92  POJ^l  JCHP(lXOJC)*YCVS(M2>P(IXD4CrYDlF(M^ 
D093J-2M2 

DO  93 1=2X2 

P(U.  l)=(P(U^)*ZCVS(3H>(U,3)*ZDIFa)yZDIF(3) 

93  P(UJ^1HP(UWZCVS(N2>WW2DIF(N1))/ZD1F(N2) 
D094K=2J^2 

P(1.1JCH>(2,1JCH>(1XK>P(2XK) 

94CXM'mNUE 

D095*=2XC 

P(U.1H>(U1WU^>P(2JX) 

950C»mNUE 
DO  96 1=2X2 

P(I.l.lHU2.IH»r[.UHU2,2) 

96CCM'mNUE 

P(l,l,lHP(l.U>+^P(lXlH>(Xl.l)y3.0 

P(LlJ4HP(L2,l,l)+P(LUl)^P(Ll.U)y3.0 

P(l,lJ^lHP(J.lJ«>+P(lXNlH>(XlXJl)y3.0 

p(iX4ia)=(P(i,M2,i>+P(ixiu)+paMi,i)y3.o 

PXlJ^14HP(UMlJ)+PXlXC,l>+P(Ll>IU)y3.0 

P(lX41JWl)^2>IlJ4l>+P(l>!2Xri>tP(l>llXI2)y3.0 

PXlMlXIl)’*aUXMlXll)H^UlX!2XllHXLlMlW 

PREF=P(IPREFJPREFja>REF) 

D097K=1XH 

D097J=1X11 

D097MX1 

97P(UJC)-P(UJ0WEF 

90CC»'mNUE 

C 

c  PRINT50 
wrile(16,50) 

IEND-0 

301  IF(1ENDXQX1)GOTO310 
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IBEO-lENIHl 
END-IENIH? 
IEha>M[NO(IENDJLl) 
c  PRINTSO 
wite(16,50) 

c  PRINT  51.(U-IBEGJEND) 
wrilB(16.5  IXlHbegiend) 
IF(MODE.EQ.3)  GO  TO  302 
c  PRINT  52.(X(D4=IBEG4END) 
wrilB(16.52Xx(i)4-a>^4cn(l) 

GO  TO  303 

c  302  PRINT  53,(X(D.I=IHEGJEND) 

302  wnte(16.53Xx(i)j=i)%4end) 

303  GO  TO  301 
310JENIM) 

c  PRINTSO 
wrile(16,30) 

311  IF(1END£QM1)  GOTO  320 
JffiG*JENI>l 
JENI^JENDf? 
JEN[>4^faN0(IENDAtl) 
c  PRINT50 
write(16,50) 

C  PRINT  54.(J4=JBEGJEND) 
wiilB(16,54XijTbcgjend) 
c  PRINT  55.(Y(J)J-jBEGJEND) 
wrile(16.55XyO')j’^b^jend() 
G0T0311 

320  KENDO 
c  PRINT50 

wnle(l6,S0) 

321  IF(KEND£Q.N1)G010330 
KffiG-KENDfl 
KEND-KENIH7 
KENDO4IN0(KENDJ41) 

c  PRINT50 
write(16.5( 

c  ?mrT  5<  _K=KBEG  JCEND) 
>vTite(16^AiUN(bqgJ(end) 
c  PRINT  57.(Z(K)4C=KffiGJKEND) 
wnte(16,S‘^z(X)4Hdb^J(end) 

GO  TO  321 
330CCKnNUE 
C 

D0999NFO,NGAM 
c  <k)999nP*5.5 

IF(.NOT.LPRINT(NF))  GO  TO  999 
c  PRINT50 
wrile(16.S0) 

c  PRINT  10,TnLE(NF) 
wnle(16,10)litle(n0 


D0998K-1J41 
c  PRINT50 
wnie(16,50) 
c  PRINT  59JC 
wTiie(16.S9)k 
1FST»1 
JFST-1 
KFST-1 

IF(NFHJ.1.0RNF£Q.4)  IFST=2 
IF(NF£Q.2.0RNF^Q.4)  JFST*2 
IF(NF£Q.3.(».NF.EQ.4)  KFST=2 
IBEO-lFST-7 
110a»mNUE 
IBE(j=IffiG+7 
IENI>=IBEG+6 
1END»MIN0(IENDJ.1) 
c  PRINT50 
wnle(16,50) 

c  PRINT  20.(1,I=IBEGJEND)  r 

wiiteC  16,20XU=ibeg,icnd) 
c  PRINT30 
write(16,30) 

JFL=jFST+Ml 

D0115JJ=JFSTM1 

J=JFL-JJ 

c  PRINT  40  J,(F(UJCJff)4“IBEGjEND) 
v/ntB(  16,40]^(fl;ijJ(.nfMNbegjend) 

115CX>KnNUE 
IF(IENDLTLl)GOTO  liO 
998CX)NTINLE 
999CXKnNUE 
RETURN 
END 
C 
C 

cxxxx:ccxxxxxxxxx(xccccxxxxxxcccxxxxmxxxmxxxxxxx(xcccccccc(xrc 

C  PROBLEM  I*PENreNTPC»TlW  C 

cxxxxxccxxxxxxccoxxxxxxxxxxcxxxxxxxcxxxxxxxxxxccccxxxxxxxxxxxrcc 

BLOCK  DATA 

c  IXXHCALLSTO* 
c  imE(BR*4NOW(14) 
c  <X)MMC»iOmASTOP4CALL4STOP 
LOaCAL  LSOL>WRIOTJLBIJaSTOP 
COMMON  F(35^.33.5XP(3545,35XRH0(35^,35XGAM(35^.35), 

1  C0N(35^^5XAKP(35^.35XAKM(35^,35XAP(35^,35). 

2  AIP05,2535)4UM(35;25.35XAJP(35^  J5),AJM(35^.35) 
CX)MMONddh(35^,33).<lelliO(35^,35Xepo(35J5,35), 

3  X(35);aJ(33XXDIF(35XXCV(33XXCVS(35),lk,c|^iaieiil. 
4Y(35).YV(35),YDIFp5),YCV(35),YCVS(35Xtiiidt,lpicv(35^,35). 

5  Z(35)ZW(35)^IF(35)W(35),ZCVS(35)^5^35), 
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6  YCM<05XYCWS(35).ARX(35).ARXK35),ARXJP(35)a^)U35a5,35X 

7  R(35)JIMN05XSX(35),SXMN(35)^CVI(35WCCV^ 

8  YCVJ(35XYCVJP05);rN«(35),ZC\rta\35)j^ 

COMMON  DU(35a5^5)X)V(35^,35XDW(33^.35XFV(35XFVP(35), 

1  FX(35)J^35)J^(35)JYM05)JT0  Wa5XTOLD(35^.35X 

2  FZ(35XravK35).VHAT(35a5.35).WHAT(35;i5^5).UOLD(35^.35) 
a:Mt04/I^roX/RELAX(13)JLJW^(13U£LJ^ 
ll^VE(10XTIMEJDT^YUl,SJ«aXW7ERO.TT^ 
2Mr>lFMAXNP>IRHO>IGAMJ.li-2aJMlM2>13>Il.N2J«. 
3ISTJST4CSTjrERJLAST, 

4IPREFJPREFjaPREFMCX£ 

OOMMON/HEADINmiLE 

CHARACTER*!©  TnUE(13) 

COMMONOmASTOPJCALUSTOP 

C(»1MC»<mA^lAEPSUfPSVJEPSWJEPSTJC(»WJim.TtX35a535)£hB^ 

1  UO(35J5.35XV0(35a5,35XW0(35;S5^5XrreRL 

CCMM0N«ESID/RMAX(13)JRES1D 

COMMON/SCX!C/SMAX.SSUM 

ooniinonface/xfoi«(35J5,35),yfoioe(35^3SX2fcice(3Sj5^S) 

C(»4MCM^/CCKF/FLOWJ)IFFACOF 

DIMe^SKW  U(35;Z5,35).V(35^  J5),W(35a5.35XPC(35^.35) 

DIMENSION  T(33^,35) 

EQUIVALENCE(F(1,1,1,IXU(1.1.1)X(F(1.1.UXV(1.1.1)X 

1  (F(l.l.l,3XW(l.l.l)X(F(l,1.1.4XPC(l.l.l)) 

2  ,(F(l,i,1.5XT(l.l.l)) 

DATA  NFMAX^IP J«HO J«5AM/5,6,7,8/ 

DATA  LSTDPi^VEJLJPRlNT/24*f ALSEy 
DATA  MODE.TIMEJIERJRESIDJTERyi.0^0.0.1/ 

DATA  RELAX.?mMESri3*.3.10*5/ 

DATAiajc/i  true; 

DATA  DTJPREF  JPREFjCPREFRHOC(»!/l£+10.1,l.l,l.(V 

(j***************************************************************** 

C  NFM^  STAND  FOR  U,V  AND  WVELOCrnES. 

C  NF=4  IS  FOR  PRESSURE 
C  NF-SISFO^lEMPERA'njRE 
C  LSOLVE-TRUE  SOLVES  THAT  PARTICULAR  PH 
DATA  (LSOLVE(IXI“l.<>y6*.1RUEy 
C  DATA(LSOLVE(IXI-1.6y4*£ALSE,.TRUE..£ALSEy 

C  IJ1UNIXNF)-TTIUE  PRINK  VARIABLE  ASSOCIATED  WriHNFCW  CALLING  PRINT 
C  DATAIJWNr(lXLPRINT(5y2*'mUEy 
C  TERMINATE  TIERATTONS  AT  TIER-LAST 
DATALAST/lOO/ 

C  UNDERELAXATTON  FACTORS 

c  DATARELAX(1)JRELAX(2XRELAX(3)£ELAX(4>0.7.07.0.7..7/ 
c  DATARELAX(5)JlELAX(6y.7,.7/ 

C  TITLES  FOR  THE  FIELD  PRINTOUTS 

DATA  TniE(lXTTILE(2XTriLE(3XTITLE(5XTTILE(6y 
lU.V.W.T.T’/ 

C  NUMKR  OF  SWEEPS  IN  THE  LINE«Y-LINETI»4AALO(»riHM 
c  DATANnMES(4XNnMES(5y2*2/ 
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C  TOLERANCES  FOR  CX3NVEIUZNCE 

DATA  EPSUJPSVJEPSW^PST/1.0E-3.1.QE.3.1.0E-3.5.0E-2/ 

DATAZEROOy 

END 


^**««**«««««««**«**«*«*««««*****«««*«******««««**«********** 

SUBROUTINE  USE 

^M********************************************************* 

C  IFUVRC£RNUMBER<^GRIDPCXNTSISTOBEUSED.THEDIMENSIC»J 
C  STATEME^^^SMUSTBEa^ANaED  THROUGH  THE  I1MXJRAM  TO  AIXOMODATE 
C  VAUJESaU\ErER1HAN(28^;»)ETC. 

Q$m000mm**0m»**m0*«*0«««m0***********»»******************» 


c  LOaCALLSTOP 
c  INIEaER*4NOW(14) 
c  (X)MMONOniASTOP,ICALL4STOP 
LCKHCAL  LSOLVEJlWNTLBlJaSTOP 
COMMON  F(35^.35,5)JV5^^5)J«0(35^^5XGAM(35^.35), 

1  CON(35^.35).AKP(35a5.35).AKM(35^,35)^5^.35). 

2  AII^5.25.35).AIM05^,35)>UP05^.35)AJM(35^.35) 
CXJMMON(Jea<35^45X<ieDiO(3545.35Xe|isi^ 

3  X05)jaJ05pa)IF05);(CV05)^CVS05Xlk.cpgalaieni. 

4  Ya5).YVa5),YDIF(35).YCV(35).YCVS(35).tinelt.lprev(35^.35). 

5  2O5)7WO5)7D1F(35)2CV(35>,ZCVS(35)j|j0(35^;}5). 

6  YCVR(35).YCVRS(35).ARX(35),ARXJ(35XARXJP(35)^1(35^.35), 

7  R(35)JRMN(35XSX(35).SXMN(35);CCVK35)^CV1P(35X 

8  YCW5),YCVIP(35).ZCVK(35XZCVKP(35)^5^^5) 

COMMON  DU(33a5^5)J)V(35^^5)DW(35^35XFV(35XFVP(35X 

1  I^5)JXN^5)JYO5)JFYM(35)JTO5XQTO5XT0LDO5^35X 

2  F2(35)JZM(35XVHAT(33^45),WHAT(35^35XU0LD(35^;J5) 
OC)MMON/INDX/REUO((13)JLFRINni3)JLaJ((ll)jn^ 
lLSOLVE(10XTTMEDT^YlJl^SRHOO»l7ERO.TLASr. 
2^ffJ4FMAXJ4PJ4RHOJ4GAMJ.lX2J-3MI^C^OJJlWO. 

3ISTJST4CSTJIERMST, 

4IPR£FJFREF,KPREFA40DE 

COMMCX>VHEADINaTIl£ 

CHARACIER*10  TnLE(13) 

OOMMWOrilASTOP4CAlIJSTa» 

COMMON«XJNVl/EPSUJEPSV£PSW,EPSr4CONVJIERl.TO(35^.35XENBAL, 

I  U0(35a535),VO(35^^5XW0(354535XnERL 

OOMMON/RESID/RMAX(13XlRESID 

Ca4MC»4/S(XlC/SMAX,SSUM 

ooinmoi^fcice^xfi»a805JS.3S).yfi)ioe(35J5.3S),zfcioe(3535.35) 

COMMONCOEF/FLOWJ)IFF^COF 

DIMENSION  U(35^.35XV(35^,35),W(35^,35XPq35^.35) 

DIMENS04  T(35^,35)jiiiaK3545.35Xxul(35X  yvl(35)^l(35) 
EQUIVALENCE(F{I.1,I,1XU(1.1,1)X(F(1,1,UXV(I,1,1)X 

1  (F(l,l,l,3XW(l,l,I)X(F(l,I.1.4)4>qi,l,l)) 

2  .(F(1.I.1,5XT(1,1.1)) 

EflMENS10NFRTIM(10) 

DATAIPR4W.TCX4IPR/1,0.38E-5,1/C*********************************************** 

ENIRYMESH 
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o  o 


DOM/UN  LENOmS  IN  THE  3  DKECnONS 

YL-5.W 
ZL-S.04 

C  NUMBBlC)F(lUDIOD^  INTHE  3  DIRECnONS 
Ll-28 
Ml-24 
Nl-24 

C 

XU(2>-0. 

XU(3)-0.08 

XU(4H).085 

XU(5H).100 

XU(6H).1274 

XU(7H).131 

XU(8>^.n4 

XU(9>4).239 

XU(10)-0J9 

XU(H)K).30 

XU(12H).33 

XU(13>^.36 

XU(14)i4).30 

XU(15)-0.65 

XU(16)K).8 

XU(17>«1.00 

XU(l8)-l.7S 

XU(19)-2.5 

XU(20)=3.25 

XU(21)-4.0 

XU(22)-4.5 

XU(23)-4.75 

XU(24)-4.8S 

XU(25)^.95 

XU(26)-5.03 

XU(27)-5.07 

XU(28)-3.10 

C 

ZW(2)K). 

ZW(3H>.48 

ZW(4>K).S8 

ZW(5H).88S 

ZW(6)-0.94 

ZW(7)^.995 

ZW(^1^ 

ZW(9)f1.396 

ZW(10W.106 

ZW(Il)-2.206 

ZW(I2)-2.511 

ZW(13)-2.566 

ZW(I4)-a.621 
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ZW(15)-2.M6 

ZW(16)-3.026 

ZW(17>-3.736 

ZW(18)-3.836 

ZW(19>-4.141 

ZW(20)-4.196 

ZW(21)-4^I 

ZW(22)-4.36l 

ZW(23)-4.66l 

ZW(24)^.04 

C 

YV(2H). 

YV(3>^.48 

YV(4>-0.58 

YV(5)^.885 

YV(6)H).94 

YV(7)-<).995 

YV(8)=1.296 

YV(9)-1.396 

YV(10)^.106 

YV(ll)-2.206 

YV(12)-«2.51l 

YV(13>-2.566 

YV(14)-2.621 

YV(15)«2.926 

YV(16)-3.026 

YV(17)»3.736 

YV(18)-3.836 

YV(19)«4.141 

YV(20)-4.196 

YV(21)»4.25l 

YV(22)-4.361 

YV(23)M.661 

YV(24)-5.04 

C 

B^YBECHN 
CINmALUME 
c  wiOTE(**)'iNrnALmiE' 

TIME-0.0 

CINmALUME  STEP 
c  WRnE(**yiNrnAL  TIME  STEP 
DT-30. 

C  riERATlONS  STOP  ATOER-LAST 
c  WRnE(*,*)inNISHnME' 

TLAST^1000.0 

C  HOW  MANY  TIMES  SHOULD  FULL  DATA  BE  PRINIED  TO  FILE 
c  WIUTE{**)’NUNfflEROFTIMESFaiPIUNnNGDATATOFlLE 
NNPIW 

C  WHENISDATATOBEFRINIED 
c  WRnE(**)’PRTIM14TmM2,...' 

PRTIM(1H).0 


PRITMO^O 

cPRTIMp^lOOO 

C  READ  RAYLEIGH  NUMBER 
c  WRnE(*.*yRAlJ-’ 

RA-1.15e3 

C  READ  FRNDIL  NUMBER 
c  WRnE(*.*)TOANTL«' 

FR-24 

C  READ  XLHC  RATIO 

c  VVRrrE(**)X3IIP  HEIGHT  TO  PACKAGE  LENCjIH  RATIO 
XLHCK)043 

c  WRnEr.*yaffl»  LENGTH  TO  PACKAGE  LENCrm  RATIO 
XLL00.017 

C  READ  RATIO  OF  CXMXXnrvrnES 
c  VVRnE(**)TUTTOOFCONDUCnTNTT^(CHIP-TX>FLUID)f 
RC-23«).0 

c  VWOTE(**)lUTTO<»^Ca®UCTTVmES(SUBSTRATE-TC>FLUID)‘ 
RS-333.0 

c  WRnEr.*)TU^TTOOFOONDUCnvrriESa‘ACKAGE-TOFLUn))* 
RI^^.0 

c  WR^m^•)TUVTK)OFCX)NDUC^^V^I^ES(UD-TOfLUIDy 
RLP271.0 

c  WRnE(*.*)lUT100F00NDUCnvrnES(AIR*TOfLUIDy 
RR-0.42 

c  WRnE(*,*yiUT100FCX)NDUCn\TnES(S0LDER-TOFLUIDy 
RO-T96.0 

c  WRrmrWTTOOFCONDUCnvniES(QOLDCOATING-TX>FLUIDy 
RK^3900.0 

C  READ  RATIOOFTHERMAL  INERTIA 
c  WRIIE(*.*)1UT100FRH0*CP(SUBSI1lATE-TOFLUlDy 
RHOCS^I.63 

c  WRnE(**)lU^TX)OFRHOC^(CHIP-TO^LUlDy 
RHOCXM).9 

c  WRriE(*,*yRATTOOFRHO^<PACKAGE.TOFLUIDy 
RHOa^l.68 

c  WRnE(*.*)lUTTOOFRHO<<y(LII>TX>fUM>y 
RHOCl^l.98 

c  WRnE(*,*yRATX)OFRHOn3»(AIR-TOFLUIDy 
RHOCRf«.00064 

c  WRnE(*,*yRATTOOFRHO^(SOLDER-TOFLUlDy 
RH0CO^.6t 

c  VW<nBr.*)lU^TTOOFRHO^P(OOLD‘TX>FLUlDy 
RHOCM-1.36 

CREAD  WHETHERTO  RAMP  SOURCE  TERMS 
c  WRnE(**yRAMP  SOURCE  TERMS  (O/iy 
IRAMI^l 

C  DETBIMINE  IF  WANT  TO  USE  A  PREVIOUSLY  OOMPIHED  SOLUTION  AS 
C  AN  INTIIAL  GUESS 

c  \MinEr.*)TIEADFROMINPUrF1LE((Viy 
IREAD-l 
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IF(IREAD£Q.  l)LREAI>.im;E 
1F(IREADJBQ.0)LREAI>-PALSE. 

C  READ  IN  RELAXATION  PARAMEHRS 
c  WRriE(*«)’ENIERRlBa«ian' 

RELAX(l)-0.4 

RELAXaH>4 

RELAX(3M).4 

RELAX(4)i«.4 

RELAX(SH).4 

RELAX(6H).4 

C  PROVIDE  INTDAL  GUESS.  TNE  FROOtAM  SOL  VES  FOR  THE  INIERIOR  peons 
CONLY.  HENCE  THE  BOUNDARY  VALUSEPQRTHETEMPBUTURES  AT  THE  HOT 
C  AND  COLD  BOUNDARIES  HAVE  BEEN  ALREADY  SPECIFIED. 

DO100I-1J.I 

DOIOOJ-IMI 

DO100K-1J41 

1MAT(UJC)^ 

RHO(UJCp>1.0 
100  CONTINUE 
DO  101 1-ai 
DOIOIMMI 
DOIOIK-IJ^I 
U(UX)-0. 

VOJJO^. 

W(UX)-0. 

T(UJCp4). 

lF(LNELl)T(UJC)-0.075 
C  SET  UP  MATERIAL  TYPE  ARRAY 
C  SUBSRATE 
]F(LLE.2)THEN 
IMAT(UJC>-.1 
RHO(UX)-RHOCS 
ENDIF 

C  COMPONENT  DETAILS 
IF  (LLE8AND.LGE4 AND. 

1  (J.GE4ANDJLE9.0RJ.GE11ANDJ1E16. 

2  ORJ.GE.18.ANDjri£.23)AND.( 

3  ICGE4AND.K1E.9.QRKGE11.ANDXLE16. 

4  ORJC.GE18.ANDJXED))THEN 
C  CHIP  AND  AIR  GAP 
IF(LEQ.7)THEN 

IF( 

1  (J.GE4ANDJIE.9.0RJ.GE11.ANDJ1E16. 

2  0RJ.GE18AM3J1E23)AND.( 

3  K.(X.4AND.IU£.9.0RJCGE11.ANDXLE16. 

4  (3RJC.GE.18.ANDJUE23))THEN 

IMAT(U4C)^3 

RHO(UJC)-RHOCR 
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ENDIF 
C  CHIP 
IR 

1  (J.GE.6.ANDJiE7.C3iU.CE.l3J^NDJlEI4. 

2  OiU.GE.20.ANDJlE21).AND.( 

3  KC£6.AI®.Ki£.7.0IUtCE13.ANDXlJE.l4. 

4  ORJC.GE.20.ANDJi-E21))™EN 
IMATXUJCy-l 
RH0(UJK)-RH0CC 

ENIXF 

EhDIF 

c  AIR  SPACE  ABOVE  CHIP 
IF(LEQ.8)THEN 
IF( 

1  (J.GE4ANDJiJE.9.0RJ.GEU.ANDJl£.16. 

2  OiU.C£.lS.ANDJl£J3)JkND.( 

3  K.GE4.AND.  KJLE9.0RJCGE  1  lANDJCLE  16. 

4  ORXGE 1*. AND. J1E23))  THEN 
IMAT(UJC)-3 

RHO(UJC)-RHOCR 

ENDIF 

ENDIF 

C  GCXIVTUNGSIENCOATINGBELOWDIECCHIP) 
IF(LBQ.6)THEN 
IF( 

1  (J.GE4.ANDJi£.9.GRJ.GE.U.ANDJI£.16. 

2  QiU.CE.18JyhB>JlEJ3VAND.( 

3  K.GE4J^ND.ia£9.0RJCGEIl.ANDi:UE16. 

4  OaK.C£.18JVNDJUE^))THEN 
IMAT(UX)-7 
RHO(UJOHIHOCM 

ENDIF 
ENDIF 
C  PACKAGE 
IF(IMAT(U4O£Q0)  THEN 
IMAT(UJC)-2 
RHCXUJC)HIHOCP 
ENDIF 
ENDIF 
C  UD 

IF(LBQ.9)THEN 

IF< 

1  (I.GE4.ANDJI£.9.CliU.GEllANDJl£.16. 

2  QIU.GE.18JkNDJI£23XAND.( 

3  ICGE4J^.IU£.9.0RXGE11ANDJCLE16. 

4  ORXCE18JU4DJi£J3))THEN 
IMAT(UJC)-4 
RHO(UJC>-RHOCL 

ENDIF 

ENMF 

C  FLUID  BEWTEEN  CHIP  AND  SUBSTRATE 
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lF(lBQ.3)riHEN 

m 

1  J.CE.4.ANDJI£.9.0IU.GE.H.ANDJLE.16. 

2  OiU.GE.l&ANDJlE23.AND. 

3  ICC£4.AND.KiJE.9.0IUCGEll.ANDJa£.16. 

4  OiUCGE.18.ANDJI£23)THEN 
IMAT(UJC)-6 

RHCXUJC)-1.0 

ENDIF 

ENDIF 

c  SOLDER  CXSNNECnCMBEWIlEN  CHIP  AND  ^JESnUlB 
IFaEQ.3)lhEn 
IF(J.EQ.4.AND. 

2  (KB).S.0iUCJQ.8.ORXH).12.OiUCEQ.lS. 

2  0iUCEQ.19.0iUCEQ.22))lben 
IMATOJJO-S 
RHO(UJKHtHOCG 
ENDF 

IF(JEQ.9.AND. 

2  (K^.5.0RXEQ.8.GRJCEQ.12.0RXEQ.15. 

2  OIUCEQ.19.CiUCEQ.22)yniEN 
IMAT(UJC)-5 
RHO(UJC)-RHOCG 
ENDIF 

IF(JEQ.11.AND. 

2  (ICEQ.5.0IUCHJ.8.0IUCEQ.12.0RJCEQ.I5. 

2  OIUCEQ.19.C»UCEQ.22)yiHEN 
IMATOJXH 
RHO(UX)-RHOCG 
ENIHF 

IF(J£Q.16.AND. 

2  aCEQ.5.0IUCGQ.8.0iUCEQ.12.0RJCEQ.lS. 

2  0RJCEQ.19.0RKEQ^)nHEN 
IMATCUJQ-4 
RHO(UJC)-«HOCG 
ENDIF 

IF(J£Q.18AND. 

2  (K^.5.0IULEQ.8.0IUC£Q.110RJCEQ.1S. 

2  0IUCEQ.19.0RJCEQ^))n«EN 
IMAT(UJO-5 
RHO(UJO-RHOCG 
ENDIF 

IFd-EQ^-AND. 

2  (KJBQ.5.(RK^.8.0iUCH).12.0RXEQ.lS. 

2  0iUCBQ.19.0RXEQ^))lNEN 
IMAlXUJO-5 
RHCKUJCHIHOCG 
ENDIF 

IFa.EQ.5.AND. 

2  (ICEQ.4.0iUCEQ.9.0RiL£Q.11.0RJCEQ.I6. 

2  0iUCEQ.I8.0IUCEQ^)yniEN 
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lMAT(UJC)-3 

RHO(UJCHU«XX} 

ENDIF 

IF(J£Q.8.AND. 

2  (K^.4.0iUCEQ.9.0iUC£Q.llQRJC£Q.16. 
2  CXUaEQ-lSORJCEQ^^lHEN 
IMAT(UJC)-5 
RHO(UJC)-RHOCG 
ENDIF 

IF(J^.12.AND. 

2  (lCEQ.4.aRJCJ3Q.9.C)fUC£Q.ll.at^^ 

2  0iUCEQ.18.0RJCBQ.23))nHEN 
lMAT(UX)-5 
RHO(UJC)-RHOCG 
ENIMF 

1F(JBQ.15.AND. 

2  (ICEQ.4.0iUC^.9.0iUC^.110iUCEQl^^ 

2  0iUCEQ.18.0IUCEQ.23))'niEN 

IMAT(UJCH 

RHCXUJCHIHOCG 

ENDIF 

IF(JflQ.19.AND. 

2  (ICEQ.4.0iUCH).9.QiUC£Q.11.0RJCEQ.i6. 
2  0RJCEQ.I8.0RJCEQJ3)nHEN 
IMAT(UJC)-5 
RHO(UJC)-RHOCG 
ENDIF 

IF(]£Q12JkND. 

2  (K^.4.0iUC£Q.9.0IUCJQ.11.0iUCJ^ 

2  0IUC£Q.18.0RJCBQ.23))>niEN 
IMATOJJC)-^ 

RHOCUJCHWOCG 

ENCHF 

ENDIF 

101  CONTINUE 
doS76i-Ul 
doS76j>‘l4nI 
dDS76lE*l^I 
wrte(20,*)  yJwlio(ijJt) 

S76oQntinuB 

IFCNQTIREAD)  RETURN 
C  READ  DATA  FROM  iraVT  FILE 
REWIND  (7) 

READ(7,*)  XYZXUl,YVl^l.U.V,Wf ,T 
ICHK*« 

DO102MXI 

lF(ABS(XUl(D-XU(I)).(jr.  lE-7)  THEN 
ICHK-1 

WRTIEI*  •)  1XU(IXXU1(D’J.XU(D^I(D 
ENDIF 


110 


102  COKITNUE 
DO103J-2M1 

IF(ABS(YVl(J).YV(J)).Gr.  lE-7)  THEN 
ICHK»1 

WRnE(*.*)  'J.YV(J).YV1(J)’J.YV(J).YV1(J) 

ENrap 

103  CX)NmNUE 
DO104K-2>Il 

IF(ABS(ZW1(K>YV(K)).CT.  lE-7)  THEN 
ICHK>1 

WRTTEC*  *)  KYV(K)ZWl(KyjC,YV(K)ZWl(K) 

ENDIF 

104  CCM^TINUE 
FaCMCEQ.DSTOP 
RETURN 

C 

ENTRY  VARRHO 

^««**«««««««««*««««««««*«««««««*««*«««»««««««*« 

RETURN 

C 

C  INGORPCIRATE  BOUNDARY  CXM)rTK»lS 
ENTRY  BNDRY 
D0864I»1X1 
D0864J-1M1 

C  ADIABATIC  SIDES  (Z-0;J=1) 

T(U.1>-T(UJ) 

T(U^1)-T(U^) 

864  CONTINUE 
D086SJ-1M1 
D0865K-1J<1 

C  ADIABATIC  SIDE  (X»0) 

T(UJQ-T(2JJC) 

T(LUJC)-0. 

865  CCNITNUE 
D0866I-1;L1 
D0866K-1.N1 

C  ADIABATIC  TOP  AND  BOTTOM  (Y"0,Y=XL) 
TaiX^TOAK) 

TOMIXH'OXCX) 

866  CONTINUE 
RETURN 

C 

ENTRY  PRTOUT 
IF(TTERNE0)GOTO400 
PRINT  401 

401  FORMAT  (IXSIMPLERV/) 

ANUCLI>0.0 
400  CONTINUE 

COMPUTE  AVERAGE  NUSSELT  NUMBER 
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ANUHaW).0 

ANUCIi>0.0 

ACHPIM). 

ACHPS-0. 

ACHPST=^. 

ACHPSBK). 

AOffB-O. 

ASUB^. 

AGAIM). 

KrCAiM).0 

HTCAPl-OO 

C  CX»mUBimWS  TO  HOT  WAliNUSSELTm  FROM  SIDES  OF  CMP 
DO  665 1-3.9 

IFaGT.3  AND.  IJLT.9)  RRT-RP 
1F(LEQ.9)RRT-RL 

C  OCKITUBimC»IS  FROM  Y-SITKS  OF  CHIP 
D0663K-11,18 
IFaEQ.3)'rHEN 
IFOMATa  10JC)JBQ.5)  THEN 

AGAI^AGAP+-2.*a(L10JC)-Ta9JC))*XCV(irZCV(Ky 
I  ((YCV(10yRG^YCV(9))) 

ELSE 

AGAP-AGAPf2.*CraiOJC>T(I.9JC))*XCV(I^ZCV(Ky 
I  ((YCV(10)+YCV(9))) 

ENDEF 

IF(IMATa.l7X).ElQ.5)  THEN 
AGAP-AGAIH^2.*(T(L17JC).T(ll8JC))*XCVa)*2CV(Ky 
1  ((VCV(17yRGfYCV(18))) 

ELSE 

AGAP^A(JAPf2.*Cr(L17JC)-T(ll8JC))*XCV(D*ZCV(Ky 
1  ((YCV(17>fYCV(18))) 

ENDIF 

ELSE 

ACHPSB-ACHPSB+-2.*(r(L10jC>.Ta9X))*XCV(D*ZCV(Ky 
1  ((YCV(10yRRT+YCV(9))) 

ACHPST-ACHPST+2.*a'ai7JC)-Tai8JC))*XCV(I)*ZC?/(Ky 
I  ((YCV(17yRRT+YCV(18))) 
e<03IF 

663  CX)NnNUE 

C  (X}NnUBini()NSFOn^Z-SIDESOFCHIP 
DO 664  >10,17 
1F(LEQ.3)T1SN 
IF(IMAT(U,  1 1)£Q.S)  THEN 

aCjAp  -agap  +2.*cr(u.i  i>t(u,  io))*Ycv(j)*xcv(iy 
1  ((ZCV(iiyRozcv(io))) 

ELSE 

AGAP  -AGAP  +2.*a’(U.10>.T(U.9))^CV(J)*XCV(iy 
I  ((ZCV(ll>fZCV(10))) 

ENDIF 
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IF(IMAT(U.  18)£Q.5)  THEN 
AGAP-AGAP4-  2.*(T(U.  18>T(U,  19))*YCV(J)*XCV(iy 
1  ((ZCV(18yRG+ZCV(19))) 

ELSE 

AGAP“AGAIH  2.*Cr(U18>-T(U.19))*YCV(J)*XCV(iy 
1  ((ZCV(18>fZCV(19))) 

ENDIF 

ELSE 

ACHPS  =ACHPS  +2.*(T(U.  1 1>T(U,  10))*YCV(J)*XCV(iy 

1  ((ZCV(liyRRT+ZCV(10))) 

2  +2.*(T(U,l8>-T(U.19))*YC:V(J)*XC:V(iy 

1  ((ZCV(18VRRT+ZCV(19))) 

ENDIF 

664  CXWnNUE 

665  CXM'mNUE 

ANUHOT-ANUHOTMOIPS+ACHPST+ACHPSB+AGAP 

D0666J=2M2 

D0667K=2J^2 

DO  669 1=2  AC 

IF(IMAT(UJC).NE.0)  THEN 

mx:APl=HTCAPl4«H0(UjC)*XCV(I)*YCV(J)*ZCV(lC)* 

1  (T(UJC>-T()(LJX)yDT»(PR'*RAr*0.5 
ENDIF 

HTCAP=inCAPHWO(UX)*XCV(I)*YCV(J)^CV(^^^ 

1  (T(U4C)-TO(UJC)yDT*(PR*RAr*0.5 
669  CONTINUE 

COMPUIE  AVERAGE  NUSSELT  NUMBER  AT  THE  HOT  AND  COLD  WALL 
IF(IMAT(4JX)£Q.0)'raEN 

C  CONIRIBtmONTOHOTWALLNUSSELTNUMBERFROMSliBSTRA'rEHEATLOSS 
ANUHar»ANUHtDT+2.*aUJJO-T(3JjC)rYCV(J)*ZCV(Ky 
1  ((xcvayRs+xcv(3))) 

•  ASUB  =ASUB  +2.*(Tl(2JX)-T(3JX))*YCV(J)*ZICV(Ky 
1  ((XCV(2yRS+XCV(3))) 

ELSE 

C  a»411UBUTK:»4PC»MTOPOFCHIP 
ANUHOT-ANUHOrr+2*(r(9JX>'T(10JJO)*YCV(J)*ZCV(Ky 
1  ((XCV(9yRL+XCV(10))) 

ACHPT  -ACHPT  +2.*(T(9  J  KVT(10JX))^CV(J)*ZCV(ICy 
1  ((XCV(9yRL+XCV(10),, 

C 

IF(IMAT(3  JX)H5.5)  THEN 
RRTi«G^C 
RRTI^GfRS 
ELSE 
RRT-l/RC 
RRTl-l/RS 
ENDIF 
C 

ACHFB  -ACHPB  +2.*RC*(r(4JJKyTa  JX))*YCV(J)*ZCV(Ky 
1  ((XCVl^HXCVOyRRT)) 

ASUBB  »ASUBB  +2.*RS*Cr(3  JX>-T(2JJOrYCV(J)*ZCV(Ky 
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1  ((XCV(3yRRTl+XCV(2))) 

ENDIF 

ANlKlI>ANUaj>+<T(L2JJt>T(LlJJC))*YCV(J)*ZCV(Ky 
1  (XDIF(L1)) 

667  CWnNUE 
666  OOKHNUE 
c  print*  •ANUHOrr.ANUHOT 
C  ANUHOT-ANUHOT 

C  MONTItXlSSlMSMAX  AND  OTHER  QUA^^TIES  AS  riERATK^jSFIlC)^^ 
C  (R>I(XH4VERGENCE,SMAXSHOULX>  BE  VERY  SMALL  (LESS  THAN  l.(£4)4) 

C  SSUM  SHOD  ACHIEVE  A  SMALL  VALUE  WTIHIN  A  FEW  riERATlC^IS.  WELL 
C  BEFC»ECXR4VERGENCE.  SSUM  WILL  NOT  BE  SMALL  IF  THE  BOUNDARY 
C  CCM)mt»«  ARE  NOT  WRITTEN  CORRECILY. 

XnER-FLOATCrrER) 

c  WRnE(21*)nER.T(22,15,l5),T(23.15.15XSSUM 
XnER5>fljOAT(rrER/IPR) 

IF(ABS(XnERm)AT(IPR)-XnER5)LT.  lE-S)  THEN 
IRESII>>1 

c  PRINT  403 JIERJTERL,'nME.SSUM.T(6.14,l4XV(6^14XENBAL 
c  WRriE(124121)TTME,ANUHC)T.ANUCLD>\SUB4nx:AP,irrCAPl 
PDIF 

405  FORMAT  (I6J4,5E12.4) 

404  FC»MAT(I6,4E12.4) 

GOviPUTE  LOCAL  AND  GLCfiAL  NUSSELT  NUMBERS  AS  riERATK^  PROCEED 
IF(  nEREQ.250.C».nEREQ.500.C».nEREQ.7S0 
LOR.riEREQ.  lOOO.ORJIEREQ. 1250.0RnERJEQ.  1500 
2.0R.riERJQ.1750.0RriEREQ.2000.C)R.nERJEQ^ 
3.ORnER^:W00.ORITERJElQ.2750.OR.riERBQ.3000 
4.0RnERBQ.3250.0RriEREQ.3500.0R.nEREQ.3730 
S.CXLTIEREQ.4000.ORIIERBQ.42S0.ORJIEREQ.4500 
6.0RriEREQ.4750.0R.riEREQ.5000.C».riEREQ.5250 
7.0RriEREQ.5500.0RnEREQ.3750.0RJrEREQ.6000 
8.C».nEREQ.62SO.ORriEREQ.6SOO.CXLriER£Q.67SO 
9miTERB3.7000.aLISTW.GT.0.ORTTMECE.TLAST)THEN 
WRriE(8.*)TIERRMAX(l)RMAXa)TlMAX(3)RMAX(Sy 
WRnE(8.*)nERRMAX(l)RMAX(2)RMAX(3)RMAX(S) 

WRnE(8*)’  V(10,14,14XT(8.14,14XV(6^14XT(10.14,14)’ 

WRnE(8,*)  V(10,14,14XT(8.l4,14XV(6a2,14XT(6A14) 
WRnE(8,*)'ANUHOT=>^,ANUHOT,*  HTCAPl-’JTTCAPl 
WRnE(8,*)fANUCLD='>U4UCLD,'  HTCAP-’JfTCAP 
WRnE(8.*)'ACHPS-'^CHPS,'  ACHPTACHPT 
WRnE(8,*)’ACHPB“>CHPB;  ASUB’^^SUB 
WRnE(8,*)'ACHPSB“'>^CHPSB,'  ACHPST’^CHPST 
WRnE(8.*)’ASUro-'.ASUBB,’  AGAP’^GAP 
C 

IFCnME.GETLAST.«LLSTOP)  THEN 
WRnE(8.*yRA-'RA,'  PR*' 

WRriE(8.*)'ZL='^L, '  RCyRC 
WRrrE(8  •yRS-'jis,’  xlho'xjk 
ENDOF 
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ENDIF 


C 

C  CHECK  TO  SEE  WHEraER  TO  WRire  DATA  TO  FILE 
IF(TTMEGEJPRTI\I(IIPR))  THEN 
WRTIEC?  *)  X.YZXU.YV;ZW.U.  V.  WP.T 
npR-npR+i 

WRriE(8,*)  •  DATA  WRITTEN  TO  FILE  FC»  TIME  =='.TIME 
ENDIF 
C 

C  SEE  WHETHER  TIMESTEP  SHOULD  E£  LENGTHENED 
C 

IF(rrERLLE.3.AND.nER.GT.  1)  THEN 
i9tiine.LE.10.)llien 
DT-DT*1.5 

WRnE(8  *)  TIME-  '.TIME,'  TIME  STEP  CHANCED  TO  DT  =•  J)T 
dse 

DT=€)T*I.l 

WRITECS  *)  TIME-  '.TIME,'  TIME  STEP  CHANCED  TO  DT  -'XT 
endif 
ENDIF 
C 

r(TlME.CE.TLAST.C».ISTC»*.CjT.O)  THEN 
C  CET  A  FIELD  PRINTOUT  AFTER  ITERATIONS  STOP 
CALL PRINT 

C  WRITEINC»IERTOOCX^4MENCENEWPRC»LEM 

WRIIE(7.*)  X,YZXU,YV.ZW,U,V,WJ>,T 

WRIIE(8.*) '  DATA  WRITTEN  TO  FILE  FOR  TIME -'.TIME 

C  WRIIEMATWXCECONDUCTTVITIESTOAFILEFDRUSEBYPIjOTROUTINE 

DO  11221-1X1 

D01122J-1X11 

D01122K-1J«II 

IF(1MAT(UK)£Q  -1)  THEN 

PC(LJK>-RS 

ELSEIF(IMAT(UK).EQ.  1}  THEN 
PC(UJO-RC 

ELSEIF(IMAT(U2C)EQ.2)  THEN 
PC(UJCH«* 

ELSEIF(IMAT(UK).EQ-3)  THEN 
FQUJCHIR 

ELSEIF(IMAT(UK)-EQ-4)  THEN 
PCOJJCHIL 

ELSEIF(1MAT(UJC)EQ.S)  THEN 
PC(UJO-RG 

elseif(imat(UJO-eq.7)  then 

PC(UKHIM 

H-SP 

PCftJJCM 

ENDIF 

1122  CCJNTTNUE 
REW1ND(9) 

WRnE(9,*)PC 
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ENDIF 

1121PC3RMAT(F13.1.5G13.5) 

RETURN 

Q^m***m*m*m*******m****m0***m*m************m»*** 

ENTRY  DIFFUS 

IF  (NFHJ.4)  RETURN 

PR-24.0 

RA-1.15e3 

RC-2360.0 

R02360.0 

RS-333.0 

RP-266.0 

RL-271.0 

RR><).42 

RO-796.0 

RM-3900.0 

IX)  500 1-1X1 

DOSOOJ-lMl 

D0300K-1^1 

C  DIFFUSIVriYPORTHEU.VORWEQUAT10NS 
IF(IMAT(UX).NE.O)  lUEN 

C  SETDIFFUSIVnYTOAHICHVALUEPORALLSOLIDREGIONS 
GAM(UX)-10E1S 
ntia? 

GAM(UJCHPR«A)**0.5 

ENIMF 

C  DTFFUSTVrTYFORTOEENERGYBQUATION 
IF(NFEQJ)THEN 
IF(IMAT(UX)XQ  *1)  THEN 
GAM(UJC)-RS/<RA*PR)«0.5 
ELSE1F(IMAT(UJ0-EQ-  1}  THEN 
GAM(UJC)-Ra(PR*RAr*03 
ELSEIF(IMAT(UJC)£Q.2)THEN 
GAM(UJC)HTP/pi*RA)**0.3 
ELSEIF(IMATXUJC)£Q.3)  THEN 
GAM(UJC)-RR/(PR*RAr*0.5 
a5ElF(IMAT(UJQXQ.4)  THEN 
GAM(UX)-RL/(PR*RAr*0.5 
ELSEIF(IMAT(UJC)EQ.S)  THEN 
GAM(UJC)«ROI(PR*RA)«0.5 
ELSEIF(IMAT(U JQEQ.7)  THEN 
GAM(I4XH»W*RAr*T).5 
FIAF. 

GAM(UX>-lAPR*RAr*0.5 

ENDIF 

C  SPECIFY  ZERO  DlFFUSIVniES  FOR  THE  ADIABATIC  BOUNDARIES 
GAMaiJQ^O 
GAMOMIJCHI.O 
GAM(UJC)-0.0 
GAM(UJ41)-0.0 


116 


GAM(U.1H).0 

ENDIF 

500  CXKTINUE 

IF  (NF.NE2  AND.  NF.NE5)  RETURN 
C  SOURCE  TERMS  ARE  EVALUATED  CM. Y  FOR  imERKXlPaNrS 
DOSOltOO^ 

DO501J-2M2 

DOS01K-2>Q 

IF(NFEQ.2)THEN 

C 

C  SOURCE  TERM  FOR  X  MOMENTUM 
C 

C  INTERPOLATE  TO  GET  THE  VALUE  OF  TEMPERATURE  AT  THE  CONIRCt 
C  VOLUME 

C  INIERFACE,  SINCE  THE  VELOCriYU  IS  EVALUATED  AT  THE  INTERFACE. 

TM-FY(J)*TXUJCHT(M(J)*T(U-lJfO 

CONl-TM 

C  SUPPLY  CM-Y  PART  <y  THE  SOURCE  TERM  IN  THE  MOMENTUM  EQ.  TO 
C  AVOID 

C  DIVEROENCE  BEFORE  200 ITERATKKSIS 
CON(UJCH=LOAT(rTER**2)*CONl/200.**2 

IF  (nERGT.200.ORIRAMPElQ.0  KXMWJKHXWl 
C 

C  S0URCETERMF0RENERGYBQUAT10NINSLAB 
C 

ELSE 

IF(IMAT(UJC)B3. 1)THEN 
0C»4a4JK>“V(PR*RA)**0.5«LHOXLLC**2 
ENDIF 
ENDIF 

501  CONTINUE 
RETURN 

ENTRY  NRGBAL 
C 

ANUOJMl.O 

HTCAP-OO 

D07»)J=2M2 

DO780K-2J« 

HTCAP''>«TCAFHRHO(UJC)*XCV(I)*YCV(J)*ZC^ 

1  (T(UJC>-TO(U4C)yDT*(PR*RA)**0.5 
770  CONTINUE 

ANUCLD=ANUai>+<TaJJJO-T(LUJC))*YCV(J)*ZCV(Ky 
1  (XD1F(L1)) 

780  CONTINUE 

C  WR^TE(^*)’H^CAP-■^^X:AP.'ANUCLD='U^NUCLD 
cprinl*>VNUCLD4nCAP 
ENBAL-ANUOXHHTCAP-l 
RETURN 
END 
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(DWriNLIE 

c  ww^E(*,•)•HrcAP-’i^x::Ap;ANlXlD-‘.A^^^^ 

c  iJrinl*,ANUQ-DHTCAP 
ENBAL-ANUaimnCAP-l 
RETURN 
END 
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